{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 数值积分"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Newton-Cotes求积公式\n",
    "\n",
    "- trapezium & simpson"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trapezium_integral(f, a, b):\n",
    "    \"\"\"梯形求积公式\n",
    "    \"\"\"\n",
    "    return (b - a) * (f(a) + f(b)) / 2\n",
    "\n",
    "\n",
    "def simpson_integral(f, a, b):\n",
    "    \"\"\"辛普森求积公式\n",
    "    \"\"\"\n",
    "    return (b - a) * (f(a) + 4 * f((a + b) / 2) + f(b)) / 6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Newton-Cotes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def costes_coefficient(n, k):\n",
    "    \"\"\"求柯特斯系数\n",
    "    \"\"\"\n",
    "    ckn = ((-1) ** (n - k))\n",
    "    ckn /= n * factorial(k) * factorial(n - k)\n",
    "\n",
    "    h = 1\n",
    "    t = Symbol('t')\n",
    "    for j in range(n+1):\n",
    "        if j != k:\n",
    "            h *= (t - j)\n",
    "\n",
    "    ckn *= integrate(h, (t, 0, n))\n",
    "\n",
    "    return ckn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "打出一张表测试科特斯系数："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1/2\t 1/2\t \n",
      "1/6\t 2/3\t 1/6\t \n",
      "1/8\t 3/8\t 3/8\t 1/8\t \n",
      "7/90\t 16/45\t 2/15\t 16/45\t 7/90\t \n",
      "19/288\t 25/96\t 25/144\t 25/144\t 25/96\t 19/288\t \n",
      "41/840\t 9/35\t 9/280\t 34/105\t 9/280\t 9/35\t 41/840\t \n"
     ]
    }
   ],
   "source": [
    "for i in range(6):\n",
    "    for j in range(i+2):\n",
    "        print(costes_coefficient(i+1, j), end= \"\\t \")\n",
    "    print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def newton_cotes_integral(f, a, b, n):\n",
    "    \"\"\"牛顿-科特斯求积公式\n",
    "    \"\"\"\n",
    "    step = (b - a) / n\n",
    "    xs = [a + i * step for i in range(n+1)]\n",
    "    return (b - a) * sum([\n",
    "        costes_coefficient(n, k) * f(xs[k]) \n",
    "        for k in range(0, n+1)\n",
    "    ])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "牛顿-科特斯求积公式测试："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "r = newton_cotes_integral(sin, 0, pi, 10)\n",
    "\n",
    "x = Symbol('x')\n",
    "print(abs(r - integrate(sin(x), (x, 0, pi))) < 0.00001)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 复化求积公式\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 复化梯形公式"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def composite_trapezium_integral(f, a, b, epsilon, max_iter=10000):\n",
    "    \"\"\" 复化梯形公式\n",
    "\n",
    "    Args:\n",
    "        f: 要求积的函数\n",
    "        a, b: 求积区间\n",
    "        epsilon: 目标精度，达到则停止，返回积分值\n",
    "        max_iter: 最大迭代次数，超出这个次数迭代不到目标精度，则 raise 一个 RuntimeError\n",
    "\n",
    "    Returns:\n",
    "        i, iter\n",
    "        i: 最终得到的积分值\n",
    "        iter_times: 迭代次数\n",
    "\n",
    "    Raises:\n",
    "        RuntimeError: 无法在 max_iter 步内迭代到目标精度\n",
    "    \"\"\"\n",
    "    m = 1\n",
    "    h = b - a\n",
    "    t = h * (f(a) + f(b)) / 2\n",
    "\n",
    "    _iter_times = 0\n",
    "    \n",
    "    t_next = 0\n",
    "    for _iter_times in range(int(max_iter)):\n",
    "        h /= 2\n",
    "        s = sum([f(a + (2 * k - 1) * h) for k in range(1, m+1)])\n",
    "        t_next = t / 2 + h * s\n",
    "        m <<= 1\n",
    "        if abs(t_next - t) <= epsilon:\n",
    "            break\n",
    "        t = t_next\n",
    "    else:\n",
    "        raise RuntimeError('无法在 max_iter 步内迭代到目标精度')\n",
    "\n",
    "    return t_next, _iter_times+1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "测试复化梯形公式："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "result_composite_trapezium_integral = (0.333343505859375, 7), actual: 1/3\n"
     ]
    }
   ],
   "source": [
    "def f(x):\n",
    "    return x ** 2\n",
    "\n",
    "x = Symbol('x')\n",
    "a = integrate(f(x), (x, 0, 1))\n",
    "\n",
    "t = composite_trapezium_integral(f, 0, 1, 0.0001)\n",
    "print(f\"result_composite_trapezium_integral = {t}, actual: {a}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 复化 Simpson 公式"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def composite_simpson_integral(f, a, b, epsilon, max_iter=1e6):\n",
    "    \"\"\"复化 Simpson 公式\n",
    "\n",
    "    Args:\n",
    "        f: 要求积的函数\n",
    "        a, b: 求积区间\n",
    "        epsilon: 目标精度，达到则停止，返回积分值\n",
    "        max_iter: 最大迭代次数，超出这个次数迭代不到目标精度，则 raise 一个 RuntimeError\n",
    "\n",
    "    Returns:\n",
    "        i, iter\n",
    "        i: 最终得到的积分值\n",
    "        iter_times: 迭代次数\n",
    "\n",
    "    Raises:\n",
    "        RuntimeError: 无法在 max_iter 步内迭代到目标精度\n",
    "    \"\"\"\n",
    "    m = 1\n",
    "    h = (b - a) / 2\n",
    "    i = h * (f(a) + 4 * f((a+b) / 2) + f(b)) / 3\n",
    "    \n",
    "    _iter_times = 0\n",
    "\n",
    "    i_next = 0\n",
    "    for _iter_times in range(int(max_iter)):\n",
    "        h /= 2\n",
    "\n",
    "        s0 = sum([f(a + (2 * k - 1) * h) for k in range(1, 2 * m + 1)])\n",
    "        s1 = sum([f(a + (4 * k - 2) * h) for k in range(1, m + 1)])\n",
    "        i_next = i / 2 + h * (4 * s0 - 2 * s1) / 3\n",
    "        \n",
    "        m <<= 1\n",
    "        if abs(i_next - i) <= epsilon:\n",
    "            break\n",
    "        i = i_next\n",
    "    else:\n",
    "        raise RuntimeError('无法在 max_iter 步内迭代到目标精度')\n",
    "\n",
    "    return i_next, _iter_times+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "result_composite_simpson_integral = 0.0416666666666667*sin(1) + 0.174637027764799, actual: 3/4 - cos(1)\n",
      "diff < 0.0001: True\n"
     ]
    }
   ],
   "source": [
    "def f(x):\n",
    "    return 3 * x ** 3 - 2 * x + sin(x)\n",
    "\n",
    "x = Symbol('x')\n",
    "a = integrate(f(x), (x, 0, 1))\n",
    "\n",
    "t, times = composite_simpson_integral(f, 0, 1, 0.0001)\n",
    "print(f\"result_composite_simpson_integral = {t}, actual: {a}\\ndiff < 0.0001: {abs(t-a) < 0.0001}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "龙贝格算法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def romberg_integral(f, a, b, epsilon, max_iter=1e6):\n",
    "    \"\"\"Romberg 积分\n",
    "\n",
    "    Args:\n",
    "        f: 要求积的函数\n",
    "        a, b: 求积区间\n",
    "        epsilon: 目标精度，达到则停止，返回积分值\n",
    "        max_iter: 最大迭代次数，超出这个次数迭代不到目标精度，则 raise 一个 RuntimeError\n",
    "\n",
    "    Returns:\n",
    "        result, T, iter\n",
    "        \n",
    "        result: 最终得到的积分值\n",
    "        T: 计算过程表\n",
    "        iter_times: 迭代次数\n",
    "\n",
    "    Raises:\n",
    "        RuntimeError: 无法在 max_iter 步内迭代到目标精度\n",
    "    \"\"\"\n",
    "    m = 0\n",
    "\n",
    "    T = [[None] * 4]\n",
    "    T[0][0] = (b - a) * (f(a) + f(b)) / 2\n",
    "\n",
    "    for m in range(1, int(max_iter)):\n",
    "        T.append([None] * 4)\n",
    "        \n",
    "        h = (b - a) / (2 ** m)\n",
    "        \n",
    "        _s = sum((f(a + (2 * k - 1) * h) for k in range(1, 1 + 2**(m-1))))\n",
    "        T[m][0] = T[m-1][0] / 2 + h * _s\n",
    "        \n",
    "        _t = 4      # 4 ** 1\n",
    "        T[m-1][1] = (_t * T[m][0] - T[m-1][0]) / (_t - 1)\n",
    "\n",
    "        if m > 1:\n",
    "            _t *= 4  # 4 ** 2\n",
    "            T[m-2][2] = (_t * T[m-1][1] - T[m-2][1]) / (_t - 1)\n",
    "        if m > 2:\n",
    "            _t *= 4  # 4 ** 3\n",
    "            T[m-3][3] = (_t * T[m-2][2] - T[m-3][2]) / (_t - 1)\n",
    "        \n",
    "        if (m > 3) and (abs(T[m-3][3] - T[m-4][3]) < epsilon):\n",
    "            break\n",
    "        # End for\n",
    "    else:\n",
    "        raise RuntimeError('无法在 max_iter 步内迭代到目标精度')\n",
    "\n",
    "    return T[m-3][3], T, m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "result_romberg_integral = 31*sin(10)/1620 + 31*sin(2)/1620 + 7392.4159348371, actual: cos(2) - cos(10) + 7392\n",
      "diff < 0.0001: True\n",
      "7 1216\n",
      "[[4*sin(10) + 4*sin(2) + 12000, 4*sin(10)/3 + 4*sin(2)/3 + 7390.50978400961, 28*sin(10)/45 + 28*sin(2)/45 + 7392.36344871677, 124*sin(10)/405 + 124*sin(2)/405 + 7392.30668995305], [2*sin(10) + 2*sin(2) + 8542.8823380072, 2*sin(10)/3 + 2*sin(2)/3 + 7392.24759467258, 14*sin(10)/45 + 14*sin(2)/45 + 7392.30757680873, 62*sin(10)/405 + 62*sin(2)/405 + 7392.36701663887], [sin(10) + sin(2) + 7679.90628050623, sin(10)/3 + sin(2)/3 + 7392.30382792522, 7*sin(10)/45 + 7*sin(2)/45 + 7392.36608789152, 31*sin(10)/405 + 31*sin(2)/405 + 7392.39496529595], [sin(10)/2 + sin(2)/2 + 7464.20444107047, sin(10)/6 + sin(2)/6 + 7392.36219664363, 7*sin(10)/90 + 7*sin(2)/90 + 7392.39451408651, 31*sin(10)/810 + 31*sin(2)/810 + 7392.40894498176], [sin(10)/4 + sin(2)/4 + 7410.32275775034, sin(10)/12 + sin(2)/12 + 7392.39249424633, 7*sin(10)/180 + 7*sin(2)/180 + 7392.40871949902, 31*sin(10)/1620 + 31*sin(2)/1620 + 7392.4159348371], [sin(10)/8 + sin(2)/8 + 7396.87506012233, sin(10)/24 + sin(2)/24 + 7392.40770542073, 7*sin(10)/360 + 7*sin(2)/360 + 7392.41582209745, None], [sin(10)/16 + sin(2)/16 + 7393.52454409613, sin(10)/48 + sin(2)/48 + 7392.41531480515, None, None], [sin(10)/32 + sin(2)/32 + 7392.6926221279, None, None, None]]\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "def f(x):\n",
    "    return 3 * x ** 3 - 2 * x + sin(x)\n",
    "\n",
    "x = Symbol('x')\n",
    "a = integrate(f(x), (x, 2, 10))\n",
    "\n",
    "t, T, iter_times = romberg_integral(f, 2, 10, 1e-10)\n",
    "print(f\"result_romberg_integral = {t}, actual: {a}\\ndiff < 0.0001: {abs(t-a) < 0.0001}\")\n",
    "print(iter_times, len(str(T)))\n",
    "print(T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def romberg_integral_sw(f, a, b, epsilon, max_iter=1e6):\n",
    "    \"\"\"Romberg 积分\n",
    "\n",
    "    用滑动窗口改写了 romberg_integral。\n",
    "    在计算的过程中，不断删除不再需要的值，\n",
    "    对于需要多次迭代的问题，可以有效节省内存。\n",
    "\n",
    "    Args:\n",
    "        f: 要求积的函数\n",
    "        a, b: 求积区间\n",
    "        epsilon: 目标精度，达到则停止，返回积分值\n",
    "        max_iter: 最大迭代次数，超出这个次数迭代不到目标精度，则 raise 一个 RuntimeError\n",
    "\n",
    "    Returns:\n",
    "        result, T, iter\n",
    "\n",
    "        result: 最终得到的积分值\n",
    "        T: 计算过程表（只有最后的 m-4 到 m 步）\n",
    "        iter_times: 迭代次数\n",
    "\n",
    "    Raises:\n",
    "        RuntimeError: 无法在 max_iter 步内迭代到目标精度\n",
    "    \"\"\"\n",
    "    m = 0\n",
    "\n",
    "    T = [[None] * 4]\n",
    "    T[0][0] = (b - a) * (f(a) + f(b)) / 2\n",
    "\n",
    "    for m in range(1, int(max_iter)):\n",
    "        T.append([None] * 4)\n",
    "        # 下面的索引 T[-1] 即 T_m, T[-2] 即 T_{m-1}, ...\n",
    "        \n",
    "        h = (b - a) / (2 ** m)\n",
    "        \n",
    "        _s = sum((f(a + (2 * k - 1) * h) for k in range(1, 1 + 2**(m-1))))\n",
    "        T[-1][0] = T[-2][0] / 2 + h * _s\n",
    "        \n",
    "        _t = 4      # 4 ** 1\n",
    "        T[-2][1] = (_t * T[-1][0] - T[-2][0]) / (_t - 1)\n",
    "\n",
    "        if m > 1:\n",
    "            _t *= 4  # 4 ** 2\n",
    "            T[-3][2] = (_t * T[-2][1] - T[-3][1]) / (_t - 1)\n",
    "        if m > 2:\n",
    "            _t *= 4  # 4 ** 3\n",
    "            T[-4][3] = (_t * T[-3][2] - T[-4][2]) / (_t - 1)\n",
    "        \n",
    "        if (m > 3) and (abs(T[-4][3] - T[-5][3]) < epsilon):\n",
    "            break\n",
    "\n",
    "        if (m > 4):\n",
    "            T = T[1:] # 清除不需要的\n",
    "        # End for\n",
    "    else:\n",
    "        raise RuntimeError('无法在 max_iter 步内迭代到目标精度')\n",
    "\n",
    "    return T[-4][3], T, m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "result_romberg_integral_sw = 31*sin(10)/1620 + 31*sin(2)/1620 + 7392.4159348371, actual: cos(2) - cos(10) + 7392\n",
      "diff < 0.0001: True\n",
      "7 850\n",
      "[[sin(10) + sin(2) + 7679.90628050623, sin(10)/3 + sin(2)/3 + 7392.30382792522, 7*sin(10)/45 + 7*sin(2)/45 + 7392.36608789152, 31*sin(10)/405 + 31*sin(2)/405 + 7392.39496529595], [sin(10)/2 + sin(2)/2 + 7464.20444107047, sin(10)/6 + sin(2)/6 + 7392.36219664363, 7*sin(10)/90 + 7*sin(2)/90 + 7392.39451408651, 31*sin(10)/810 + 31*sin(2)/810 + 7392.40894498176], [sin(10)/4 + sin(2)/4 + 7410.32275775034, sin(10)/12 + sin(2)/12 + 7392.39249424633, 7*sin(10)/180 + 7*sin(2)/180 + 7392.40871949902, 31*sin(10)/1620 + 31*sin(2)/1620 + 7392.4159348371], [sin(10)/8 + sin(2)/8 + 7396.87506012233, sin(10)/24 + sin(2)/24 + 7392.40770542073, 7*sin(10)/360 + 7*sin(2)/360 + 7392.41582209745, None], [sin(10)/16 + sin(2)/16 + 7393.52454409613, sin(10)/48 + sin(2)/48 + 7392.41531480515, None, None], [sin(10)/32 + sin(2)/32 + 7392.6926221279, None, None, None]]\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "def f(x):\n",
    "    return 3 * x ** 3 - 2 * x + sin(x)\n",
    "\n",
    "x = Symbol('x')\n",
    "a = integrate(f(x), (x, 2, 10))\n",
    "\n",
    "t, T, iter_times = romberg_integral_sw(f, 2, 10, 1e-10)\n",
    "print(f\"result_romberg_integral_sw = {t}, actual: {a}\\ndiff < 0.0001: {abs(t-a) < 0.0001}\")\n",
    "print(iter_times, len(str(T)))\n",
    "print(T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "题目 4.4 实验内容和步骤\n",
    "\n",
    "$$\n",
    "I_1 = \\int_0^1 e^{-x^2} dx\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "integral exp(-x ** 2) dx from 0 to 1\n",
      "\n",
      "actual result (by sympy): sqrt(pi)*erf(1)/2 = 0.746824132812427\n",
      "\n",
      "result_trapezium=0.6839397205857212\n",
      "result_simpson=0.7471804289095104\n",
      "result_composite_trapezium=0.7468238989209475, times_composite_trapezium=9\n",
      "result_composite_simpson=0.7468241406069852, times_composite_simpson=4\n",
      "result_romberg=0.7468241326473878, times_romberg=4\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "def f1(x):\n",
    "    if not isinstance(x, float):\n",
    "        x = float(x)\n",
    "    return math.exp(-x ** 2)\n",
    "\n",
    "result_trapezium = trapezium_integral(f1, 0, 1)\n",
    "result_simpson = simpson_integral(f1, 0, 1)\n",
    "result_composite_simpson, times_composite_simpson = composite_simpson_integral(f1, 0, 1, 1e-6)\n",
    "result_composite_trapezium, times_composite_trapezium = composite_trapezium_integral(f1, 0, 1, 1e-6)\n",
    "result_romberg, _, times_romberg = romberg_integral(f1, 0, 1, 1e-6)\n",
    "\n",
    "actual = integrate(exp(-Symbol('x') ** 2), (Symbol('x'), 0, 1))\n",
    "\n",
    "print(f'''integral exp(-x ** 2) dx from 0 to 1\n",
    "\n",
    "actual result (by sympy): {actual} = {actual.evalf()}\n",
    "\n",
    "result_trapezium={result_trapezium}\n",
    "result_simpson={result_simpson}\n",
    "result_composite_trapezium={result_composite_trapezium}, times_composite_trapezium={times_composite_trapezium}\n",
    "result_composite_simpson={result_composite_simpson}, times_composite_simpson={times_composite_simpson}\n",
    "result_romberg={result_romberg}, times_romberg={times_romberg}\n",
    "''')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "I_2 = \\int_0^1 \\frac{\\sin(x)}{x} dx\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "integral sin(x) / x dx from 0 to 1\n",
      "\n",
      "actual result (by sympy): Si(1) = 0.946083070367183\n",
      "\n",
      "result_trapezium=0.9207354924039483\n",
      "result_simpson=0.9461458822735868\n",
      "result_composite_trapezium=0.9460829746282349, times_composite_trapezium=9\n",
      "result_composite_simpson=0.9460830853849476, times_composite_simpson=3\n",
      "result_romberg=0.9460830703672595, times_romberg=4\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def f2(x):\n",
    "    if not isinstance(x, float):\n",
    "        x = float(x)\n",
    "    return math.sin(x) / x\n",
    "\n",
    "result_trapezium = trapezium_integral(f2, 1e-32, 1)\n",
    "result_simpson = simpson_integral(f2, 1e-32, 1)\n",
    "result_composite_simpson, times_composite_simpson = composite_simpson_integral(f2, 1e-32, 1, 1e-6)\n",
    "result_composite_trapezium, times_composite_trapezium = composite_trapezium_integral(f2, 1e-32, 1, 1e-6)\n",
    "result_romberg, _, times_romberg = romberg_integral(f2, 1e-32, 1, 1e-6)\n",
    "\n",
    "actual = integrate(sin(Symbol('x')) / Symbol('x'), (Symbol('x'), 0, 1))\n",
    "\n",
    "print(f'''integral sin(x) / x dx from 0 to 1\n",
    "\n",
    "actual result (by sympy): {actual} = {actual.evalf()}\n",
    "\n",
    "result_trapezium={result_trapezium}\n",
    "result_simpson={result_simpson}\n",
    "result_composite_trapezium={result_composite_trapezium}, times_composite_trapezium={times_composite_trapezium}\n",
    "result_composite_simpson={result_composite_simpson}, times_composite_simpson={times_composite_simpson}\n",
    "result_romberg={result_romberg}, times_romberg={times_romberg}\n",
    "''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def adaptsim(f, a, b, eps=1e-8, max_iter=10000):\n",
    "    \"\"\"自适应 Simpson 求积\n",
    "    \n",
    "    P.S. 这个函数名来自 Gander, W. and W. Gautschi, “Adaptive\n",
    "         Quadrature – Revisited,” BIT, Vol. 40, 2000, pp. 84-101.\n",
    "         该文档可以在 https://people.inf.ethz.ch/gander/ 找到。\n",
    "         但该函数的实现并没有使用此文中的递归方法。\n",
    "\n",
    "    Args:\n",
    "        f: 要求积的函数\n",
    "        a, b: 求积区间\n",
    "        eps: 目标精度，达到则停止，返回积分值\n",
    "        max_iter: 最大迭代次数，超出这个次数迭代不到目标精度，则 raise 一个 Exception\n",
    "\n",
    "    Returns: (I, m, p)\n",
    "        I: 积分的近似值\n",
    "        m: 分层数\n",
    "        p: 分点\n",
    "\n",
    "    Raises:\n",
    "        Exception: 无法在 max_iter 步内迭代到目标精度\n",
    "    \"\"\"\n",
    "    p = [a, b]  # 分点\n",
    "    p0 = p\n",
    "    ep = [eps]\n",
    "    m = 0\n",
    "    q = 0\n",
    "    I = 0\n",
    "    \n",
    "    for _iter_times in range(int(max_iter)):\n",
    "        n1 = len(ep)\n",
    "        n = len(p0)\n",
    "        \n",
    "        if n <= 1:\n",
    "            break\n",
    "        \n",
    "        h = p0[1] - p0[0]\n",
    "        s0 = h /  6 * ( f(p0[0]) + 4 * f(p0[0] + h/2) + f(p0[0] + h  ) )\n",
    "        s1 = h / 12 * ( f(p0[0]) + 4 * f(p0[0] + h/4) + f(p0[0] + h/2) )\n",
    "        s2 = h / 12 * ( f(p0[0] + h/2) + 4 * f(p0[0] + 3*h/4) + f(p0[0] + h) )\n",
    "        \n",
    "        if abs(s0 - s1 - s2) <= 15 * ep[0]:\n",
    "            I += s1 + s2\n",
    "            p0 = p0[1:]\n",
    "            \n",
    "            if n1 >= 2:\n",
    "                ep = ep[1:]\n",
    "            \n",
    "            q += 1\n",
    "        else:\n",
    "            m += 1\n",
    "            p0 = [p0[0], p0[0] + h/2] + p0[1:]\n",
    "            \n",
    "            if n1 == 1:\n",
    "                ep = [ep[0]/2, ep[0]/2]\n",
    "            else:\n",
    "                ep = [ep[0]/2, ep[1]/2] + ep[1:]\n",
    "            \n",
    "            if q == 0:\n",
    "                p = p0\n",
    "            else:\n",
    "                p = p[:q] + p0\n",
    "        \n",
    "    else:\n",
    "        raise Exception('无法在 max_iter 步内迭代到目标精度')\n",
    "\n",
    "    return I, m, p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VPW9//HXh0UWRaGCVggQFFB2kLCJIrgi9oKKLWBcQDFqBVqtfbjkqkgb/dXbiuJSBGRR07q0FemtvVyxYBFFCYILIMoFEgKoIYgLYc/n98cJIZAJGZJJTpJ5Px+PeZA5851zPhzCZ858l88xd0dERGq+WmEHICIilUMJX0QkTijhi4jECSV8EZE4oYQvIhInlPBFROKEEr6ISJxQwhcRiRNK+CIicaJOWAdu2rSpJyYmhnV4EZFqafny5dvcvVlZ3htawk9MTCQjIyOsw4uIVEtmllnW96pLR0QkTijhi4jECSV8EZE4EVofvohUrH379pGdnc3u3bvDDkXKoH79+iQkJFC3bt2Y7VMJX6SGys7OplGjRiQmJmJmYYcjx8Ddyc3NJTs7mzZt2sRsv6V26ZjZTDP72sw+LeF1M7MpZrbOzD42s7OjOfDy5ZCYCOnpxxixiERl9+7dnHzyyUr21ZCZcfLJJ8f821k0ffizgcFHef0yoF3BIwX4Y7QHz8yElJSSk356evChUKvW0T8com0nEm+U7Kuvivi3K7VLx93/bWaJR2kyDHjeg3slLjWzxmZ2mrtvjSaAvDxIHf8dydvnQN26waNhQ9I/PIuUJ7uSt6c2cPDDwQEjOfnQ+9PTgw+NvDyKtAt+LtpORCTexWKWTgtgU5Hn2QXbijGzFDPLMLPDVlxlfXMCTJgAt90GY8fCNdeQ+vsmhcn+oLw8I/X6TdClC1x8MVx3Hak//6Yw2R9qB6mpxY+vbwIilWfKlCl06NCB5BCuvK6++mrWr19fpvcevM/3xIkTC59H2lbUypUr6devH506daJr1668/PLLha+NHDmSL774okyxxFqlDtq6+zRgGoBZUuEZa9US+DAH9u0LHjt3ktWpNUS4v3pWfgK0bQtffgmLF5P13UkRj5WVmQ/DroTTT4czziB9yyBSHu9I3q7ga5K+CYhUrGeeeYYFCxaQkJBQqcddtWoVBw4c4PTTTy/T+1NTU+nTpw+5ublMmDCBG2+8kVdeeaXYtu7duxe+p2HDhjz//PO0a9eOLVu20LNnTy699FIaN27MbbfdxqOPPsr06dNj9Vcsu4OfXkd7AInApyW89iwwqsjztcBppe+zp4N7w4buL77oxbRu7Q7FH61bH9GuVX7kdg2+cu/cOTgAeGs2RG536i73nTuLHf/FF4NjmQV/RopRpCpbvXp1aMe+5ZZbvG7dut65c2d/7LHHonrPCy+84L169fJu3bp5SkqK79+/3zdu3Oht27b1nJwcP3DggJ977rk+f/5837Bhg5955pl+zTXX+FlnneXDhw/3nQX/j++9916fNWtWxGM89thjPmbMGHd3//jjj71Tp06F7yvq1ltv9caNG/tnn3121G0l6dq1q3/++efu7n7gwAFPTEz0ffv2RXUeior0bwhkeBR5O9IjFlf484BxZvYS0Af41qPsv2/dGtLSIl9hp6Ud3jcP0LBhsP2wdg9b5HbTToHkT4K8vnUrWQmnRf7G8NVx0KgudOoEvXpBr16kb7uElEfakJenbwNSQ/zyl7ByZWz32b07PP54xJemTp3K//zP/7Bw4UKaNm3KwoULueOOO4q1a9iwIe+++y5r1qzh5ZdfZsmSJdStW5ef//znpKenc/3113P33Xdz22230bt3bzp27Mgll1zCxo0bWbt2Lc899xz9+/fnxhtv5JlnnuGuu+5iyZIljBo1KmJcv/jFLxg4cCCvvfYaaWlpPPvsszRs2PCwNv/5n//J4MGDqVOnDk8//TQ33XQTr776arFt3bp1i3iMDz74gL1793LGGWcAUKtWLdq2bctHH31Ez549j+UMx1ypCd/M/gwMBJqaWTbwIFAXwN2nAm8AQ4B1QB4wJpoD9+wJR6uddjCxpqZCVha0ahX5w6HUdmbQvDmtWgWJ+0itmu2GW1Nh2TJ4/XWYOZNUNpDH4SPkB8cFlPBFjt2gQYNYeZQPnLfeeovly5fTq1cvAHbt2sUpp5wCwNixY3n11VeZOnXqYfto2bIl/fv3B+Daa69lypQp3HXXXWzdupVmzSIXk6xVqxazZ8+ma9eu3HLLLYXvL+o3v/kNZsaKFSuYOHEi7k7Xrl2LbYtk69atXHfddcyZM4datQ4NkZ5yyimFXT1himaWTuSPykOvO3B7zCIqIjk5ugQbTbsSvzFMbgjJk4IN7rBhA1ltSxg/yMyHO++CAQPg/POhSRPS00v/UBIJXQlX4pWltCt8d+eGG27gkUceKdYmLy+P7OxsAH744QcaNWoEFJ+2ePB5gwYNCuevv/baazz00EMAzJgxg6SkJL744gtOOOEEtmzZEjHWg/s5OEBb9DiRth303Xffcfnll5OWlkbfvn0Pe2337t00aNAg4vEqVVn7gsr76Nmz5zH3Z5VXtP3yJY4f1NvqXr9+8KRWLX+x7QPesO6ew9qUNCYhUtnC7MN3d2/durXn5ORE1XbVqlXetm1b/+qrr9zdPTc31zdu3Oju7uPGjfO0tDR/8cUX/fLLL3d39w0bNjjg7777rru733TTTf773//e3d1HjBjhb775ZsTj7Nixw9u3b+9r1671iy++2F999dVy/R0P2rNnj19wwQU+efLkiK937tzZt27desz7jXUfflwl/Gi9+GLhWG/xRL57t/vixe4PPOCtj9sS+YMhYX/YfwWRapXw3d1feukl79atm3fp0sXPPvtsf++993zRokXep08f378/+D915ZVX+syZMwsHbZOTk/2ss87yq666qnDw9fnnn/fU1NSIxxgzZow/8cQT7u6elZXlZ5xxRuGHTHm88MILXqdOHe/WrVvhY8WKFe7u/uWXX3qvXr3KtF8l/EoSzbcBs8jfBIwD7v37u//+9+7r1h3TPkViJeyEX5E2bNjgnTp1ivhaXl7eYR8SYXvsscd8xowZZXpvrBO+yiOXIDkZNm6E/Pzgz0j98q1aRX5vq5O+g5074a67gjUDXbuSfuWrpIw9QGZm8LFQWlkJESmbBg0a8NBDD7F58+awQwGgcePG3HDDDWGHAagefrmkpQUDv0U1bAhpTzeGFStg/XqYPBmaNCF1bm/ydh+5cjjyimARObrExEQ+/TRiPUcALr30UlqVdEVWycaMGUOdOlWjMLESfjkkJ8O0acF6ArPgz2nTinwbaNMmmP/89ttkWeRfvqxMD6b4iIhUMCX8coqm6wegVavIle9akRl8UgwYADNnwvffq+aPiFQIJfxKUmL3z2MNghdzcuCmm0hvOoGUG/aor19EYk4Jv5KU2P1zx6lw332wejW8+y6px/2evAP1Dnuv+vpFJBaU8CvRUbt/zKBfP7J2nhzxvVmZ+fDWW8Flv0gNMXv2bMaNG1em927cuJE//elPhc8zMjKYMGFCrEI7Juecc06pbR5//HHyjqzlXsmU8KuYEqd61toMF10EHTrAE0/Ajh2AavxL7FS336UjE35SUhJTpkwJJZZ333231DZK+FJMiX39z50Kzz8PTZoEM39atCD9ghmkjM1Xf7+U28E7x8Xyd+mKK66gZ8+edOrUiWnTphVunzVrFu3bt6d3794sWbKkcPvf//53+vTpQ48ePbjooov46quvgKB+zXXXXUe/fv1o165dYV35e+65h8WLF9O9e3cmT57MokWL+MlPfkJ+fj6JiYnsKLgoAmjXrh1fffUVOTk5DB8+nF69etGrV6/Djn/Q7NmzGTZsGAMHDqRdu3aFtXgAHnvsMTp37kznzp15vEh9ohNOOAGARYsWMXDgQK6++mrOOusskpOTcXemTJnCli1bGDRoEIMGDeLAgQOMHj2azp0706VLFyZPnlz2E30syrpiq7yPqr7SNkylrshdvtz9ppu8tWVGdc8AiU/HstI22vtPHIvc3Fx3D1a+durUybdt2+Zbtmzxli1b+tdff+179uzxc845x2+//XZ3d9++fbvn5+e7u/v06dP9zjvvdHf3Bx980Lt27ep5eXmek5PjCQkJvnnzZl+4cGFhbR13P+z5hAkTfObMme7uvnTpUr/wwgvd3X3UqFG+ePFid3fPzMz0s846q1jcs2bN8h//+Me+bdu2wtiXLVvmGRkZ3rlzZ//hhx/8+++/944dO/qHH37o7u7HH398YQwnnniib9q0yQ8cOOB9+/YtPF7RUhMZGRl+0UUXFR7zm2++iXgOq2I9fImxUqt/nn02zJhB1szI/flZWcG9f0WiVdJSkPIsEZkyZQqvvfYaAJs2beKLL77gyy+/ZODAgYXli0eMGMHnn38OQHZ2NiNGjGDr1q3s3buXNm3aFO5r2LBhNGjQgAYNGjBo0CA++OADGjduXOKxR4wYwaRJkxgzZgwvvfQSI0aMAGDBggWsXr26sN13333HDz/8UHiFftDFF1/MyScH42lXXXUV77zzDmbGlVdeyfHHH1+4ffHixfTo0eOw9/bu3bvwLl/du3dn48aNnHvuuYe1Of3001m/fj3jx4/n8ssv55JLLonyrJaPunSqsRLn9ls2PPggfP11JUck1VWJY0dlXKy6aNEiFixYwHvvvcdHH31Ejx49CksWl2T8+PGMGzeOTz75hGefffaw9iWVQi5Jv379WLduHTk5OcydO5errroKgPz8fJYuXcrKlStZuXIlmzdvLpbsy3K8ourVOzTLrnbt2uzfv79YmyZNmvDRRx8xcOBApk6dytixY6Pef3ko4VdjEfv76x0grdvLMGlSMPdzwgTIyqp2A3JSuUocO0qL3L403377LU2aNKFhw4Z89tlnLF26FIA+ffrw9ttvk5uby759+3j11VcPe0+LFi0AmDNnzmH7e/3119m9eze5ubksWrSIXr160ahRI77//vuIxz94NX7nnXfSoUOHwqv1Sy65hCeffLKwXUk3ZXnzzTfZvn07u3btYu7cufTv35/zzjuPuXPnkpeXx86dO3nttdc477zzoj4nRePdtm0b+fn5DB8+nN/+9rd8+OGHUe+nPNSlU41FvttXbZKT74K1/wGPPgp//CPpT31DSq0ZhfP7dctGOVK0d5iL1uDBg5k6dSodOnTgzDPPLLwhyGmnncbEiRPp168fjRs3PuxG4BMnTuSnP/0pTZo04YILLmDDhg2Fr3Xt2pVBgwaxbds27r//fpo3b06zZs2oXbs23bp1Y/To0cW6VkaMGEGvXr2YPXt24bYpU6Zw++2307VrV/bv38+AAQOYOnVqsfh79+7N8OHDyc7O5tprryUpKQmA0aNH07t3byC4E9eRxzyalJQUBg8eTPPmzXn88ccZM2YM+fn5ABFv/FIRzEOa152UlOQZR7vHocRGVhaJnU8g8/sfFXupdetgPYDUTGvWrKFDhw5hh1FuEydO5IQTTuCuu+6qlOPNnj2bjIwMnnrqqUo53tFE+jc0s+XunlSW/ekKv6Zr1YqsHyK/pMFdkfiiK/w4kJgY+QburdnIxivugIkToVu3yg5LKlhNucKPZ7G+wtegbRyIOCDXwEm7ajksXAjdu8PVV0NBfXEN8NYcYV3QSflVxL+dEn4ciFi4bbqR/NfhQSf+Aw/A//5vcGeuvlO0ereGqF+/Prm5uUr61ZC7k5ubS/369WO636i6dMxsMPAEUBuY4e7/74jXWwMzgWbAduBad88+2j7VpVPFbN8Of/gDiY/cQqYXn3ytAd7qZ9++fWRnZ5c6/12qpvr165OQkEDdunUP216eLp1SE76Z1QY+By4GsoFlwCh3X12kzavAf7v7HDO7ABjj7tcdbb9K+FVTrVqOe/GBXLOgyqeIhKui+/B7A+vcfb277wVeAoYd0aYj8K+CnxdGeF2qiRJX7570LezaVcnRiEgsRZPwWwCbijzPLthW1EfAVQU/Xwk0MrPIhd2lSos4wFtrN2k7boO2bWH2bF3qi1RTsRq0vQs438xWAOcDm4EDRzYysxQzyzCzjJycnBgdWmIp4gDv8/VJfvtWSEiAMWMgKSmY3SMi1Uo0ffj9gInufmnB83sB3D3iWmAzOwH4zN0TjrZf9eFXQ/n58PLLcM89wfr7YcOC8g3t24cdmUjcqOg+/GVAOzNrY2bHASOBeUcE0NTMDu7rXoIZO1LT1KoFo0bBZ5/Bww8Ht1zs1Al+8QvIzQ07OhEpRakJ3933A+OA+cAa4BV3X2Vmk8xsaEGzgcBaM/scOBUoY409qRYaNIB774V16+Cmm+Cpp4L+/aeeggilYEWkalBpBSm/Tz8Nbrv41lvQpQtMmQIDB4YdlUiNpNIKEq7OneHNN+Gvf4XvvoNBg2DECNi0qfT3ikilUcKX2DCDq66C1auDYmzz5sGZZ8Jvfwu7d6s+j0gVoIQvsdWwYXB7xTVrYMgQuP9+0lvdQ8pN+1WfRyRk6sOXirVgAYmXdSBz/5Fr9VSfR6QsdAMUqbouuoisA5EvKrKyKjkWkTinLh2pcCXW5zlFVRxFKpMSvlS4iPV5LI+0r26EsWPhm2/CCUwkzijhS4WLWJ9nRi2Sf50QFGPr0CEo2aAbdYhUKA3aSrhWrICbb4bly+Hyy+GZZ6BV8RuwiEhAC6+k+urRA5YuhcmTYdEi6NgRnngCDhQrtioi5aSEL+GrUycozbBqFQwYEPzcty98/HHYkYnUKEr4UnW0bg3/+Af8+c/BnM2ePYNFXHv3hh2ZSI2ghC9VixmMHBmUaBg5EiZNChL/smVhRyZS7SnhS9V08snwwgvw97/D9u1BF8/dd+u+uiLloIQvVdtPfhL07d94Y3B3re7dYcmSsKMSqZaU8KXqa9wYpk8PSjDv3QvnnQd33EH6rL2qwClyDFRLR6qPiy6CTz6Bu+8m/fGvSLED5BUsIzlYgROChV4iUpwWXkm1lHjqLjK/blBsuypwSk2nhVcSd7Jyiid7UAVOkaNRwpdqqaTqC61O3AH79lVuMCLVhBK+VEsRK3DW3kPatz+Hfv3gs8/CCUykClPCl2opYgXOOfVI/svwoBO/Rw946inIzw87VJEqI6qEb2aDzWytma0zs3sivN7KzBaa2Qoz+9jMhsQ+VJHDJScHuT0/P/gzORkYPjyYyTNoEIwfD4MHw+bNIUcqUjWUmvDNrDbwNHAZ0BEYZWYdj2j2n8Ar7t4DGAk8E+tARaJ22mlBTZ4//jFYpNWlS1BvXyTORXOF3xtY5+7r3X0v8BIw7Ig2DpxY8PNJwJbYhShSBmZw661Bvf127YK6PMnJsGNH2JGJhCaahN8C2FTkeXbBtqImAteaWTbwBjA+JtGJlFf79sFV/kMPBVf53bvDO++EHZVIKGI1aDsKmO3uCcAQ4AUzK7ZvM0sxswwzy8jJyYnRoUVKUacOPPBAkOhr14bzz4f779f0TYk70ST8zUDLIs8TCrYVdRPwCoC7vwfUB5oeuSN3n+buSe6e1KxZs7JFLFJWffvCypVw/fXw298GNXn+7//Cjkqk0kST8JcB7cysjZkdRzAoO++INlnAhQBm1oEg4esSXqqeRo1g1qyge2ft2qCLZ/Zs3UBd4kKpCd/d9wPjgPnAGoLZOKvMbJKZDS1o9ivgZjP7CPgzMNrDKtIjEo2f/Sy4hWLPnjBmDIwcSfr0naq+KTWaiqdJfDtwAP7rv0i/bxUpPEueH1q+27BhsLhL1TelKilP8TQlfBEg8bQ9ZH5Zr9h2Vd+UqkbVMkXKKeur4skeVH1TahYlfBGOUn2zaV7lBiJSgZTwRSih+qbtIi1nLEyYAHv2hBOYSAwp4YtQQvXNWXVI/sUp8OSTcM45sG5d2GGKlIsGbUVKM28ejB4N+/cHnwojR4YdkcQxDdqKVKShQ4MVul26wKhRcPPNkKe+fal+lPBFotGqFSxaBPfeCzNmQO/esHp12FGJHBMlfJFo1a0LDz8M8+dDTg4kJQVlGUSqCSV8kWN1ySVBF0/fvkFZhjFj1MUj1YISvkhZnHYavPlmUHZ5zpygi2fNmrCjEjkqJXyRsqpdO7ixyvz58PXXQRfPiy+GHZVIiZTwRcrr4ouDLp6kJLjuumAWz65dYUclUowSvkgsNG8Ob70F990XzOLp0wc+/zzsqEQOo4QvEit16gQ1Gv75T9iyBZKSSJ+wVDX2pcqoE3YAIjXO4MGwYgXpg2aQ8mQXDs7fycyElJTgZ9XYlzDoCl+kIrRsSeq+ieRx/GGb8/IgNTWkmCTuKeGLVJCsTRZ5u2rsS0iU8EUqSIk19k/cEdxaUaSSKeGLVJCINfZr7yHt25/DkCGwbVs4gUncUsIXqSARa+zPOY7kaYOCQmxnnw0ffBB2mBJHlPBFKlBycnAT9Pz84M/kZAsWZi1ZEszVPO88mDoVQrovhcSXqBK+mQ02s7Vmts7M7onw+mQzW1nw+NzMdsQ+VJEaJCkJPvwQLrwQbrsNbrhBBdikwpWa8M2sNvA0cBnQERhlZh2LtnH3O9y9u7t3B54E/lYRwYrUKD/6Efz3fwf1eF58Mai++cUXYUclNVg0V/i9gXXuvt7d9wIvAcOO0n4U8OdYBCdS49WqFVTc/Oc/YfPm4Mr/9dfDjkpqqGgSfgtgU5Hn2QXbijGz1kAb4F/lD00kjlx6adDF0749XHFFsDpLUzclxmI9aDsS+Iu7R/xNNbMUM8sws4ycnJwYH1qkmmvdGhYvhrFjgztrXXaZpm5KTEWT8DcDLYs8TyjYFslIjtKd4+7T3D3J3ZOaNWsWfZQi8aJ+fZg+PXj8+9/QsydkZIQdldQQ0ST8ZUA7M2tjZscRJPV5RzYys7OAJsB7sQ1RJA6NHQvvvBP8fO658Nxz4cYjNUKpCd/d9wPjgPnAGuAVd19lZpPMbGiRpiOBl9w1oVgkJpKSYPlyGDAg+AC4+WbYvTvsqKQai6oP393fcPf27n6Gu6cVbHvA3ecVaTPR3YvN0ReRcmjaNJjBc/DGKgMGkD4lVzX2pUxUD1+kqqtdOyjMk5RE+qi/k/KLBqqxL2ViYfXAJCUleYYGo0SOSWKLfWRuqVtse+vWQekGqfnMbLm7J5XlvaqlI1KNZG0tnuxBNfYlOkr4ItVIiTX2T9tXuYFItaSEL1KNRKyxTx5p39wG//hHOEFJtaGEL1KNRKyxP3knyWd9CP/xHzBpUlCLWSQCDdqK1AS7dsEtt8ALL8DQofD883DSSWFHJRVAg7Yi8a5BA5gzB6ZMgTfegN69YfXqsKOSKkYJX6SmMIPx4+Gtt2DHDujTh/RffqBFWlJIC69EapoBA2D5ctIHTiPliU5apCWFdIUvUhMlJJC67yHyOP6wzXl5Qal9iU9K+CI1VNYmi7xdi7TilhK+SA1V4iKtk3dWbiBSZSjhi9RQERdp1dpF2rab4Z57dAvFOKSEL1JDRVykNbMuybeeBL/7HQwZAtu3hx2mVCItvBKJRzNmwO23Q0ICvPYadO0adkQSJS28EpFjM3YsvP12cAetfv1IH/+e5uvHAc3DF4lXffsG8/UHTCXlqa6arx8HdIUvEs9+/GNS9zyo+fpxQglfJM5pvn78UMIXiXOarx8/lPBF4pzm68cPJXyROFfifP1bTgzm619+uebr1xBRJXwzG2xma81snZndU0Kbn5nZajNbZWZ/im2YIlKRkpNh48bgZlkbN0LyDXVg6tTgk+Bf/4JeveCTT8IOU8qp1IRvZrWBp4HLgI7AKDPreESbdsC9QH937wT8sgJiFZHKdvPNwXz9XbugXz/4y1/CjkjKIZor/N7AOndf7+57gZeAYUe0uRl42t2/AXD3r2MbpoiEpl8/yMgIVuP+9Kdw333q16+mokn4LYBNRZ5nF2wrqj3Q3syWmNlSMxscaUdmlmJmGWaWkZOTU7aIRaTyNW8OCxcGK7IeeSS4Yfo334QdlRyjWA3a1gHaAQOBUcB0M2t8ZCN3n+buSe6e1KxZsxgdWkQqRb168OyzwWPBgqBf/9NPw45KjkE0CX8z0LLI84SCbUVlA/PcfZ+7bwA+J/gAEJGaJiUFFi2CnTuD8gx//WvYEUmUokn4y4B2ZtbGzI4DRgLzjmgzl+DqHjNrStDFsz6GcYpIVXLOObB8OXTpAldfHdRhUL9+lVdqwnf3/cA4YD6wBnjF3VeZ2SQzG1rQbD6Qa2argYXAr909t6KCFpEqoHnz4Ep/7Fh4+OGgX3/HjrCjkqNQPXwRKb9nn4Xx44NVW6+9Bp07hx1RjaV6+CISrltuCWbx/PCD+vWrMCV8EYmN/v0P69dPH/Yyia1dN1WpQnQDFBGJnYJ+/fTLXiBl3k/IIyi9rJuqVA26wheR2KpXj9T1Y3VTlSpICV9EYq6km6fopirhUsIXkZgr8aYqjb7RfP0QKeGLSMxFvKlKnT2kfXc7DBkCuVqmEwYlfBGJuYg3VZldj+TpFwSLtZKSYOXKsMOMO1p4JSKV6/33Yfjw4C5aM2bANdeEHVG1ooVXIlJ99OkTzNdPSgq+Ctx5J+kvHCAxEc3Zr2Cahy8ile/UU+Gtt+BXvyJ98pek1NpHXn5tQHP2K5K6dEQkVIlNfyAz94Ri21u3Du6vK4dTl46IVFtZ24sne9Cc/YqghC8ioSpxzn7LcHofajIlfBEJVcQ5++wkre5E2HzkzfWkPJTwRSRUEefsj/+E5C//AD17wuLFYYdYYyjhi0jokpODAdr8/ODP5Cl9g/n6J54IF1wATzwBIU0wqUmU8EWkaurUCZYtg8svh1/+Eq69lvSZuzVfvxyU8EWk6jrpJPjb3yAtjfQ/QcpYJzMzuNg/OF9fST96mocvItVC4qm7yPy6QbHt8TZfX/PwRaTGy8opnuxB8/WPhRK+iFQLJc7Xb6H6+tGKKuGb2WAzW2tm68zsngivjzazHDNbWfAYG/tQRSSelThff/ed8OGH4QRVzZSa8M2sNvA0cBnQERhlZh0jNH3Z3bsXPGbEOE4RiXMR5+tP3Epy/b/BOefArFlhh1jlRXOF3xtY5+7r3X0v8BIwrGKQepe5AAAJsklEQVTDEhEprth8/QfbBlf3554LN94YTNvZvTvsMKusaBJ+C2BTkefZBduONNzMPjazv5hZy5hEJyJSmmbNYP58uPdemD4dzjsvmLMpxcRq0PbvQKK7dwXeBOZEamRmKWaWYWYZOTk5MTq0iMS92rXh4Ydh7lz4/HM4+2zS7/5Ii7SOEE3C3wwUvWJPKNhWyN1z3X1PwdMZQM9IO3L3ae6e5O5JzZo1K0u8IiIlGzYMMjJIbziWlEfbapHWEaJJ+MuAdmbWxsyOA0YC84o2MLPTijwdCqyJXYgiIsegXTtS7RHyOP6wzXl5kJoaUkxVRKm3OHT3/WY2DpgP1AZmuvsqM5sEZLj7PGCCmQ0F9gPbgdEVGLOIyFFlZUe+lo33RVoqrSAiNU5iYuRx29Y/+p6N204I5nVWUyqtICJSRMRFWrV2k7b9Frj+eti5M5zAQqaELyI1TsRFWnOOI3lSh2DktndvWBN/Q43q0hGR+LJgAVxzTTCKO306jBoVdkTHRF06IiLRuugiWLECuneHa64h/eJZJLb2uJivX+osHRGRGqdFC1i4kPQrXiXljWHkEQziHpyvD0G3UE2jK3wRiU9165K66pq4mq+vhC8icaukeflZWTXzhulK+CISt0q8qcpxX0F2duUGUwmU8EUkbkWcr3/cftK4LxjU/ec/wwmsgijhi0jcijhff2Ydkj+6G5o3hyFD4L77YP/+sEONCSV8EYlrxW6qkgyceSa8/z6MHQuPPAIXXkj6U99U+3LLmpYpIhJJgwbBwqzzzyf9xrdI+Xc98gpeqq7TN7XSVkSkFInN95K59bhi21u3Dr4VVCattBURqUBZXxZP9lD9yi0r4YuIlKLE6ZtN8yK/UEUp4YuIlCLi9E3bRVrOWLjzTti7N5zAjpESvohIKSJO35xVh+TxTWHyZOjfH9atCzvMUmnQVkSkPObOhRtvhH374Nlng9LLFUiDtiIiYbniCli5MliZm5wMN95I+szdVXLOvhK+iEh5tWoFCxfC/feTPmsPKWOdzExwPzRnvyokfXXpiIjEUOKpu8j8ukGx7bGas68uHRGRKiIrp3iyh6oxZz+qhG9mg81srZmtM7N7jtJuuJm5mZXp00dEpLorcc5+s12VG0gEpSZ8M6sNPA1cBnQERplZxwjtGgG/AN6PdZAiItVFiXP2vx4Ld98d6pz9aK7wewPr3H29u+8FXgKGRWj3G+B3wO4YxiciUq1EnLP/XC2Sbz0RHn00mLP/xRehxBZNwm8BbCryPLtgWyEzOxto6e7/iGFsIiLVUrGSy2PqwR//CH/7G6xfDz16wOzZwTSeSlTuQVszqwU8BvwqirYpZpZhZhk5OTnlPbSISPVy5ZXw0UfQqxeMGUN63ykktjxQafP1o0n4m4GWRZ4nFGw7qBHQGVhkZhuBvsC8SAO37j7N3ZPcPalZs2Zlj1pEpLpKSIAFC0gf8TopH4wlM7t2pc3XjybhLwPamVkbMzsOGAnMO/iiu3/r7k3dPdHdE4GlwFB31yR7EZFIatcmdelQ8jj+sM15eZCaWnGHLTXhu/t+YBwwH1gDvOLuq8xskpkNrbjQRERqrpLm5WdlVVy/flS3OHT3N4A3jtj2QAltB5Y/LBGRmq1Vq6Abp9h2z4In58HttwfFeGJIK21FREIQcb5+g3zSur0MEybA4MGweXPkN5eREr6ISAgiztefXovkFb8OpnAuWQJdusDLL8fsmCqeJiJSFX3xBVx3Hbz/PowaRfrAaaQ+fAKZmUm4Z1hZdhlVH76IiFSydu3gnXfgkUdIf3AtKS/VIq+c1+fq0hERqarq1IH77yf11OfI84alty+FEr6ISBWX9VW9mOxHCV9EpIorqeTysVLCFxGp4iJN4SwLJXwRkSqu6BTO8lDCFxGpBg6WXIbly8u6DyV8EZE4oYQvIhInlPBFROKEEr6ISJxQwhcRiROhFU8zs++BtaEcvOppCmwLO4gqQufiEJ2LQ3QuDjnT3RuV5Y1hFk9b6+7F7nsbj8wsQ+cioHNxiM7FIToXh5hZmcsMq0tHRCROKOGLiMSJMBP+tBCPXdXoXByic3GIzsUhOheHlPlchDZoKyIilUtdOiIicaLCE76ZDTaztWa2zszuifB6PTN7ueD1980ssaJjCksU5+JOM1ttZh+b2VtmVs7aeFVXaeeiSLvhZuZmVmNnaERzLszsZwW/G6vM7E+VHWNlieL/SCszW2hmKwr+nwwJI86KZmYzzexrM/u0hNfNzKYUnKePzezsqHbs7hX2AGoD/wecDhwHfAR0PKLNz4GpBT+PBF6uyJjCekR5LgYBDQt+vi2ez0VBu0bAv4GlQFLYcYf4e9EOWAE0KXh+Sthxh3gupgG3FfzcEdgYdtwVdC4GAGcDn5bw+hDgn4ABfYH3o9lvRV/h9wbWuft6d98LvAQMO6LNMGBOwc9/AS40szLdkb2KK/VcuPtCd88reLoUSKjkGCtLNL8XAL8BfgfsrszgKlk05+Jm4Gl3/wbA3b+u5BgrSzTnwoETC34+CdhSifFVGnf/N7D9KE2GAc97YCnQ2MxOK22/FZ3wWwCbijzPLtgWsY277we+BU6u4LjCEM25KOomgk/wmqjUc1HwFbWlu/+jMgMLQTS/F+2B9ma2xMyWmtngSouuckVzLiYC15pZNvAGML5yQqtyjjWfAOGutJUSmNm1QBJwftixhMHMagGPAaNDDqWqqEPQrTOQ4Fvfv82si7vvCDWqcIwCZrv7H8ysH/CCmXV29/ywA6sOKvoKfzPQssjzhIJtEduYWR2Cr2m5FRxXGKI5F5jZRUAqMNTd91RSbJWttHPRCOgMLDKzjQR9lPNq6MBtNL8X2cA8d9/n7huAzwk+AGqaaM7FTcArAO7+HlCfoM5OvIkqnxypohP+MqCdmbUxs+MIBmXnHdFmHnBDwc9XA//yglGJGqbUc2FmPYBnCZJ9Te2nhVLOhbt/6+5N3T3R3RMJxjOGunuZa4hUYdH8H5lLcHWPmTUl6OJZX5lBVpJozkUWcCGAmXUgSPg5lRpl1TAPuL5gtk5f4Ft331ramyq0S8fd95vZOGA+wQj8THdfZWaTgAx3nwc8R/C1bB3BIMXIiowpLFGei/8CTgBeLRi3znL3oaEFXUGiPBdxIcpzMR+4xMxWAweAX7t7jfsWHOW5+BUw3czuIBjAHV0TLxDN7M8EH/JNC8YrHgTqArj7VILxiyHAOiAPGBPVfmvguRIRkQi00lZEJE4o4YuIxAklfBGROKGELyISJ5TwRUTihBK+iEicUMIXEYkTSvgiInHi/wOkXJDWBApFHAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "\n",
    "def f(x):\n",
    "    if not isinstance(x, float):\n",
    "        x = float(x)\n",
    "    return math.exp(-x ** 2)\n",
    "\n",
    "# actual result (by sympy): sqrt(pi)*erf(1)/2 = 0.746824132812427\n",
    "\n",
    "I, m, p = adaptsim(f, 0, 1, eps=1e-9)\n",
    "\n",
    "# 画图\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "xf = np.linspace(0, 1)\n",
    "yf = [f(x) for x in xf]\n",
    "\n",
    "yp = [f(x) for x in p]\n",
    "\n",
    "plt.plot(xf, yf, 'r-', label='f=exp(-x ** 2)')\n",
    "plt.plot(p, yp, 'bo', label='adaptive points')\n",
    "plt.legend()\n",
    "plt.xlim(0, 1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD8CAYAAAB6paOMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuczXX+wPHXewYxIdduxsxoowxNZBokt6jfkBJq0XRhs9o21abaVRSp2dqyiK0YG7oo3aUSSay2KCO3kJJcBjGITUMu8/798TkzjjHMMc7M98w57+fjcR5zzvf7Pd/znq8x7/l8P5e3qCrGGGMiV5TXARhjjPGWJQJjjIlwlgiMMSbCWSIwxpgIZ4nAGGMinCUCY4yJcJYIjDEmwlkiMMaYCGeJwBhjIlw5rwMoqFatWpqQkOB1GMYYU6YsXrx4h6rWLs57Qy4RJCQkkJmZ6XUYxhhTpojIhuK+124NGWNMhLNEYIwxEc4SgTHGRLiQ6yMwxpSsgwcPkpWVxf79+70OxRRDxYoViY2NpXz58kE7pyUCYyJMVlYWVapUISEhARHxOhxzElSVnTt3kpWVRb169YJ23iJvDYnIRBHZLiLfHGe/iMgYEVkrIstF5BK/fbeKyPe+x62BBLR4MSQkwJQpAX8PxpiTsH//fmrWrGlJoAwSEWrWrBn01lwgfQSTgdQT7O8E1Pc9+gPPA4hIDWAo0BxIAYaKSPVAgtqwAfr3t2RgTEmxJFB2lcS/XZG3hlR1vogknOCQrsBL6mpeLhSRaiJyDtAOmK2quwBEZDYuobwWSGA5OTB4wB7StmRAdDSUK+ceMTHuUanS0V8rV4Zq1eCMMyCI986MMSbcBaOPoA6wye91lm/b8bYfQ0T641oTQLP87Rt3V4G//vXkI8pLCnmPmjXhrLOOPM4888jzOnXc8QGaMgUGD4aNGyEuDtLTIS3t5EM0JpKNGTOG559/nksuuYQpxWz69+vXj4EDB5KYmHjC40aPHk2NGjW45ZZbjntMr169eOyxx6hfv36h+6dOncoPP/zA4MGDixVrqAuJzmJVzQAyAESSNW97XF1g9V44fBgOHYKDB2HfPtdcyPua9/jlF9i9+9jHzz/DDz/AF1/Ajh2gemwA1apB3bruERd35Ot558HvfucShwhTprhbVjk57m15t7DAkoExJ+O5557jk08+ITY2ttjn+Pe//13kMYcOHWLixIl8/fXXJzzujjvu4KmnnmLChAmF7v/oo4+4++67ixVnmaCqRT6ABOCb4+wbD/T2e70GOAfoDYw/3nHH/6xmCqoxMaqvvKLBdeiQ6k8/qS5bpvrxx6ovv6z65JOqd96pes01qk2aqNasqerSxZFH5cqqF1+s8ZW2HbMLVOPjgxynMSVo1apVnn7+7bffruXLl9fGjRvryJEjizx+79692rlzZ01KStJGjRrp1KlTVVW1bdu2umjRIlVVPf300/Whhx7SpKQkbd68uf7000+qqjpr1iy99dZbVVX14MGDmpycrHPnzlVV1UGDBulDDz2kqqqHDx/WhIQEPXjw4DGfn5ubq0lJSZqbm3vU9pEjR2rfvn1VVXX58uXaqFEj/fXXX0/+ghRDYf+GQKYG8Pu8sEcwWgTTgQEiMhXXMbxHVbeKyCzg734dxFcBDwZywvj4ErrlEh195JbQieTkwKZNriXxww+wdi388AMbl9Uq9PCNG3LhrnsgMREaNYKLL3Z9FcaEur/8BZYuDe45mzSB0aOPu3vcuHHMnDmTuXPnUqtWLebOncu99957zHExMTF88cUXzJw5k3PPPZcPP/wQgD179hxz7K+//kqLFi1IT0/nr3/9KxMmTGDIkCF8/vnnNGvmbjeXK1eOyZMnc/311zN27FhmzpzJl19+CUBUVBTnn38+y5Ytyz8+z5IlS7j44ouP6aS95557aNeuHe+++y7p6emMHz+emJiYk7tWIaLIRCAir+E6fmuJSBZuJFB5AFUdB8wAOgNrgRygr2/fLhF5DFjkO9Vw9XUcn0izZuD5mnMxMXDBBe7hJy7B3Q4qKK7CNnjxRXd7Ks9557n/EE2buq9Nmrj+CButYcxR2rdvz9ITJKOLLrqI++67j7/97W906dKF1q1bH3NMhQoV6NKlCwDNmjVj9uzZAGzdupWGDRvmH9eoUSNuvvlmunTpwoIFC6hQoUL+vjPPPJMtW7YckwhmzpxJp06djvnMqKgoJk+eTFJSErfffjutWrU6uW88hAQyaqh3EfsVuPM4+yYCE4sXWuhJTz+6jwBczkjPOAdu3AObN8OKFbBkifsra8kSeOedIwefeSakpEDz5u6RkmItB+OtE/zlXlqKahE0aNCAr7/+mhkzZjBkyBA6dOjAI488ctSx5cuXz/+LPTo6mkOHDgFQqVKlY8bcr1ixgmrVqrF9+/ajtu/fv59KlSodE8fHH3/M22+/XWjs33//PZUrV2bLli2Bf8MhKCQ6i8uKvFtVhY8aEoiNdQ//vx5++QWWLXNJITMTvvoKPvjA7ROBCy90SaFlS2jd2r22VoOJIEW1CLZs2UKNGjW46aabqFatWkCdxHkaNmzI2rVr81+/88477Nq1i/nz59OlSxe++uorqlWrBsB3331H48aNj3r/nj17OHToEDVr1jzm3Hv27OHuu+9m/vz5DBgwgLfeeovrr78+4NhCiSWCk5SWdpJ9F1WqwOWXu0ee3bth0SL48kv3+PBDmDzZ7atd2x3bujW0aeP6G8rZP5OJXCtWrOCBBx4gKiqK8uXL8/zzzwf83k6dOnHzzTcDsGPHDgYNGsScOXOoW7cuAwYM4J577uHFF19k27ZtVKpUibPPPvuo98+ePZuOHTsWeu57772XO++8kwYNGvDCCy/Qvn172rRpw5lnnln8b9YjooUNp/RQcnKyRlxhGlXXIT1/Pnz2mfv6449uX+XKLiF07OgejRtbi8GcktWrVx913zzcdevWjaeeeuq4cwQARo0aRdWqVbntttuO2t6vXz/69etHixYtSjrMk1LYv6GILFbV5OKcz/7UDAUiUL++e+T9IGZlHUkKn34KM2a47WedBR06uKTQoYO7P2WMOa4nn3ySrVu3njARVKtWLb/l4O9kbkOVZVaPIFTFxkLv3vD887BmjRuuNHGi++U/Zw784Q9unG2jRvDAAzB3Lhw4kP/2KVPc4n1RUbaIn4lsF1xwAW3atDnhMX379qVcBN+CjdzvvKyJi4O+fd1DFVauhNmz4aOPYMwYGDHC9UdceSVTagyg/5S25Oxzed5mQBtjTsRaBGWRiOsruPde+Phj2LkTpk1zLYivvmLwv+vlJ4E8OTlutJMxxhRkLYJwULkydO3qHqpsjAYKGQOwcYPCl1/BpZe6e0bGGIO1CMKPCHFxhY8qimMjtGjhbjP95S+wcGHhi/AZYyKKJYIwlJ7uZjz7i4mB9PE14eWXITnZdUK3bAn16rmlvhcvtqRgQsrkyZMZMGBAsd67fv16Xn311fzXmZmZnq0eetlllxV5zOjRo8nxX7KglFkiCENpaZCR4QYVibivGRmQ1r8y3HST60/Yvt2tj9SoEYwa5ZJD/fowZIgbpWSMT1kcgVYwESQnJzNmzBhPYvniiy+KPMYSgSkRaWmwfj3k5rqvx4wWOuMMuOUWN6t52zb497/dQnlPPOGWuUhJgbFjITvbg+hNqMirwbFhg2swBquM7HXXXUezZs1o1KgRGRkZ+dsnTZpEgwYNSElJ4fPPP8/f/v7779O8eXOaNm1Kx44d2bZtGwDDhg3j5ptvpmXLltSvXz+/nsCgQYP47LPPaNKkCaNGjWLevHl06dKF3NxcEhIS2L17d/6569evz7Zt28jOzqZHjx5ceumlXHrppUd9fp7JkyfTtWtX2rVrR/369Xn00Ufz940cOZLGjRvTuHFjRvut4VTZV/hq3rx5tGvXjuuvv54LL7yQtLQ0VJUxY8awZcsW2rdvT/v27Tl8+DB9+vShcePGXHTRRYwaNerULnYgirt+dUk9mjVrVrwFuk1wbNmi+s9/utoMoFqunKvV8MYbqvv2eR2dCYKTqUcQH39s/Y1g1ODYuXOnqqrm5ORoo0aNdMeOHbplyxatW7eubt++XX/77Te97LLL9M4771RV1V27duXXA5gwYYIOHDhQVVWHDh2qSUlJmpOTo9nZ2RobG6ubN2/WuXPn6tVXX53/ef6v7777bp04caKqqi5cuFA7dOigqqq9e/fWzz77TFVVN2zYoBdeeOExcU+aNEnPPvts3bFjR37sixYt0szMTG3cuLHu3btXf/nlF01MTNSvv/5aVV2thLwYqlatqps2bdLDhw9rixYt8j8vPj5es7OzVVU1MzNTO3bsmP+ZP//88zFxhGI9AhNOzjkHBg50j2++cX0Kr7wC778P1au7W0u33ebWQDJhb+PGk9seqDFjxvDuu+8CsGnTJr7//nt++ukn2rVrR+3atQHo2bMn3333HQBZWVn07NmTrVu3cuDAAerVq5d/rq5du1KpUiUqVapE+/btj1pIrjA9e/Zk+PDh9O3bl6lTp9KzZ08APvnkE1atWpV/3P/+9z/27t2b/xd9niuvvDJ/Ebru3bvz3//+FxGhW7dunH766fnbP/vsM5o2bXrUe1NSUvKrsjVp0oT169dzuf86ZMB5553HunXruOuuu7j66qu56qqrAryqxWe3hszxNW4M//iH+18/ezakprrOhiZNXJ/CuHFQSJEQEz6Ot4LJqaxsMm/ePD755BMWLFjAsmXLaNq06TFLRRd01113MWDAAFasWMH48eOPOr5gwZiCrwtq2bIla9euJTs7m2nTptG9e3cAcnNzWbhwIUuXLmXp0qVs3rz5mCRQnM/zd9ppp+U/918u21/16tVZtmwZ7dq1Y9y4cfTr1y/g8xeXJQJTtOhot7bRq6/Cli1uJvOBA3DHHa4Fccst8PnnNuooDB13BFp68c+5Z88eqlevTkxMDN9++y0LFy4EoHnz5vznP/9h586dHDx4kDfffPOo99SpUweAF1988ajzvffee+zfv5+dO3cyb948Lr30UqpUqcIv/oWi/OT99T5w4EAaNmyY/9f9VVddxdixY/OPO97S2LNnz2bXrl3s27ePadOm0apVK1q3bs20adPIycnh119/5d133y20gM7x+Me7Y8cOcnNz6dGjB48//niR9ZaDIaBEICKpIrJGRNaKyKBC9seLyBwRWS4i80Qk1m/fUyKyUkRWi8gYOZn0aUJPjRpw112uxsJXX7kkMG2aWzq7SRMYPx727vU6ShMkxx2BdgpLlaSmpnLo0CEaNmzIoEGD8lf2POeccxg2bBgtW7akVatWR62uOWzYMG644QaaNWtGrVpHl4xNSkqiffv2tGjRgocffphzzz2XpKQkoqOjufjiiwvtbO3ZsyevvPJK/m0hcLerMjMzSUpKIjExkXHjxhUaf0pKCj169CApKYkePXqQnJzMJZdcQp8+fUhJSaF58+b069fvmNtCJ9K/f39SU1Np3749mzdvpl27djRp0oSbbrqJJ554IuDzFFtRnQhANPADcB5QAVgGJBY45k3gVt/zK4CXfc8vAz73nSMaWAC0O9HnWWdxGfTLL6oZGUc6mKtUUR0wQHXlSq8jM4Xwunh9MA0dOlSffvrpUvu8SZMm5XdgeynYncWBtAhSgLWquk5VDwBTga4FjkkEPvU9n+u3X4GKvgRyGq7W8baTyFOmLKhcGf74R/j6a/jiC7fURUaGm6PQrp0r13n4sNdRGmOOI5BEUAfY5Pc6y7fN3zKgu+95N6CKiNRU1QW4xLDV95ilqqsLfoCI9BeRTBHJzLZx62WXiJut/PLLrp7CP/7hBp736AHnnw8jR1rnsgmqYcOGcf/995fa5/Xp04d//etfpfZ5pSVYncX3A21FZAnQFtgMHBaR84GGQCwueVwhIsf0oKhqhqomq2py3tAxU8bVru2Wrli71rUI4uLgvvtcnYW773bbKZuzVsOBWsd+mVUS/3aBJILNQF2/17G+bflUdYuqdlfVpsBg37bduNbBQlXdq6p7gY+AlkGJ3JQN0dHQrRv85z9uPaPu3d2w0wYNmHLJCPrfdjjos1bNiVWsWJGdO3daMiiDVJWdO3dSsWLFoJ63yJrFIlIO+A7ogEsAi4AbVXWl3zG1gF2qmisi6cBhVX1ERHoCfwRSAQFmAqNV9f3jfV5E1iyOND/9BM8/T8Lj/diQW/eY3fHxblkMUzIOHjxIVlZWkWP3TWiqWLEisbGxlC9f/qjtp1KzOKDi9SLSGRiNG/kzUVXTRWQ4rpd6uohcDzyB6xyeD9ypqr+JSDTwHNDGt2+mqg480WdZIogcUVGK6rGjiUWU3FwbZWzMySjxRFCaLBFEjoQEdzuooPioTaz/+6tuwlrVqqUelzFl0akkAptZbDxT6KzV0w6TnvgKDBrkMsWjj8LPP3sSnzGRwhKB8Uyhs1ZfiCZtxYNu1nKbNjBsmNvx0EO2JLYxJcRuDZnQtny5azq8+SZUqgR/+hPcf79b48gYk89uDZnwlZQEr78OK1e6iWnPPOMK6Awc6KqsGWNOmSUCUzY0bAgvveTKaPbq5RJCvXquL2HnTq+jM6ZMs0Rgypbf/Q4mTYLVq+G66+Cpp1yn8sMPg1/5QWNM4CwRmLKpQQM3BXnFCujUCR5/3CWE9HT49VevozOmTLFEYMq2Ro3gjTdg6VK30umQIa7V8NxzrniOMaZIlghMeLj4YlcgZ8ECuPBCuPNO16/w2muQm+t1dMaENEsEJry0aAFz58KMGVClCtx4o6uvPGuWldI05jgsEZjwI+L6Db7+Gl55xXUip6a6usvHqUNrTCSzRGDCV1SUm7787bduuOmyZXDJJdC3L2zeXPT7jYkQlghM+KtQ4UgxnPvug1dfdaOOhg6FvXu9js4Yz1kiMJGjWjV4+mnXQrjmGhg+3CWEF16wmsomolkiMJGnXj2YOtWNMEpIgH794NJL4bPPvI7MGE8ElAhEJFVE1ojIWhEZVMj+eBGZIyLLRWSeiMT67YsTkY9FZLWIrBKRhOCFb8wpaNECPv/cJYUdO9xqp717w6ZNXkdmTKkqMhH4qow9C3QCEoHeIpJY4LARwEuqmgQMx1Ury/MS8LSqNgRSAFspzIQOEejZ090ueuQRNxfhggvgscdg3z6vozOmVATSIkgB1qrqOlU9AEwFuhY4JhH41Pd8bt5+X8Iop6qzAXxF7HOCErkxwRQT44rgfPstdOnikkLDhvDWWzb/wIS9QBJBHcC/rZzl2+ZvGdDd97wbUEVEagINgN0i8o6ILBGRp30tDGNCU3y8W7Ji7lw44wy44Qb4v/9jyogtJCS4EakJCW6ZI2PCRbA6i+8H2orIEqAtsBk4DJQDWvv2XwqcB/Qp+GYR6S8imSKSmW1VqEwoaNcOFi+GsWOZ8t84+j9wBhs2uMbBhg3Qv78lAxM+AkkEm4G6fq9jfdvyqeoWVe2uqk2Bwb5tu3Gth6W+20qHgGnAJQU/QFUzVDVZVZNr165dzG/FmCArVw4GDGBwzfHkcPpRu3JyYPBgj+IyJsgCSQSLgPoiUk9EKgC9gOn+B4hILRHJO9eDwES/91YTkbzf7lcAq049bGNKz8bNhd/N3LjR+g5MeCgyEfj+kh8AzAJWA2+o6koRGS4i1/oOawesEZHvgLOAdN97D+NuC80RkRWAABOC/l0YU4Li4o6zXTfCE0/AwYOlG5AxQWbF640pwpQprk8gx2+8W0ylXDIaPUNa5kBo3BgyMqBlS++CNBHPitcbU4LS0tzv+fh4N+0gPh4yJkSRtuheN+9g925o1Qr+/GfYs8frcI05adYiMOZU/fKLq5k8diycdRaMGQM9erisYUwpsRaBMV6qUgVGj4Yvv3SJ4IYboGtXyMryOjJjAmKJwJhgSU6GRYtgxAj45BNXT3nCBJuZbEKeJQJjgqlcOVfzYPlyVwSnf3/o0AF++MHryIw5LksExpSE88+HOXNg/HjIzISLLoJRo6zugQlJlgiMKSlRUa5FsGoVXHEFDBwIl1/uXhsTQiwRGFPSYmPh/ffhlVfg++/dLaOnn7bWgQkZlgiMKQ0ibkLCypXQuTP89a+uEM7333sdmTGWCIwpVWedBW+/7VoHq1bBxRe7eQe5uV5HZiKYJQJjSpt/66B9e7jnHteH8OOPXkdmIpQlAmO8cu658MEHMHEiLFniRhZlZNi8A1PqLBEY4yUR6NsXVqxwi9bdfjtcey1s2+Z1ZCaCWCIwJhTExcGsWfDMMzB7tmsdTJ9e9PuMCQJLBMaEiqgouPtu+PprqFPHrVfUr59b1M6YEmSJwJhQk5joFrB78EHXf9CkCXzxhddRmTAWUCIQkVQRWSMia0VkUCH740VkjogsF5F5IhJbYH9VEckSkX8FK3BjwlqFCvD3v8P8+W5oaevWMHQoHDrkdWQmDBWZCEQkGngW6AQkAr1FJLHAYSOAl1Q1CRgOPFFg/2PA/FMP15gIc/nlsGwZ3HwzDB8ObdvC+vVeR2XCTCAtghRgraquU9UDwFSga4FjEoFPfc/n+u8XkWa4OsYfn3q4xkSgqlVh8mR49VX45ht3q+j1172OyoSRQBJBHWCT3+ss3zZ/y4DuvufdgCoiUlNEooB/4grYG2NORe/esHSp60Po1csNO9271+uoTBgIVmfx/UBbEVkCtAU2A4eBPwMzVPWEpZpEpL+IZIpIZnZ2dpBCMiYM1avn+g0efhheesktYGelXc0pCiQRbAbq+r2O9W3Lp6pbVLW7qjYFBvu27QZaAgNEZD2uH+EWEXmy4AeoaoaqJqtqcu3atYv3nRgTKcqVc/0Fc+fCvn1uItrIkTYj2RRbIIlgEVBfROqJSAWgF3DUTBcRqeW7DQTwIDARQFXTVDVOVRNwrYaXVPWYUUfGmGJo08Z1JHfp4qqiXXstU8b/QkKCm5KQkABTpngdpCkLikwEqnoIGADMAlYDb6jqShEZLiLX+g5rB6wRke9wHcPpJRSvMcZfjRrwzjswZgxTPqpO/zui2bDBNQ42bHB1cSwZmKKIhlhzMjk5WTPtnqcxJy3hnN/Y8NNpx2yPj7cRp5FARBaranJx3mszi40JExu3HZsEADZuLOVATJljicCYMBEXd5ztZ+4v3UBMmWOJwJgwkZ4OMTFHb4uRHNK393M1kkPsNrAJHZYIjAkTaWmurk18vCtzEB/vXqf1+M3VSO7WDXbv9jpME4Kss9iYcKcKY8e6IaZ168Jbb7mJaCasWGexMeb4RFydg/nz4eBBuOwymDDBbhWZfJYIjIkULVu62sht27oJBn37upnJJuJZIjAmktSqBTNmwLBhbq2iyy6DH3/0OirjMUsExkSa6GhX5OaDD9xMs2bNXL1kE7EsERgTqTp3diuX1q0LnTq58ae5uV5HZTxgicCYSPa738GCBa7WwZAh0L077NnjdVSmlFkiMCbSxcTAK6/AM8/Ahx9CSgqsXOl1VKYUWSIwxhwZYvrpp65F0KIFTJvmdVSmlFgiMMYc0bo1LF4MDRu6mciPPmr9BhHAEoEx5mh16rjJZ7fe6oaZ9ugBv/zidVSmBFkiMMYcq2JFmDQJRo2C9993k9F++MHrqEwJCSgRiEiqiKwRkbUickypSRGJF5E5IrJcROaJSKxvexMRWSAiK337egb7GzDGlBAR+Mtf3ByDrVvh0kth9myvozIloMhEICLRwLNAJyAR6C0iiQUOG4GrR5wEDAee8G3PAW5R1UZAKjBaRKoFK3hjTCno0AEWLXK3jFJTXSvB1ikKK4G0CFKAtaq6TlUPAFOBrgWOSQQ+9T2fm7dfVb9T1e99z7cA24HawQjcGFOKzjvPzTe47joYOBD++Ec4cMDrqEyQBJII6gCb/F5n+bb5WwZ09z3vBlQRkZr+B4hIClABsBuNxpRFlSvDm2+6iWcvvABXXgk7dngdlQmCYHUW3w+0FZElQFtgM3A4b6eInAO8DPRV1WPGoolIfxHJFJHM7OzsIIVkjAm6qCh47DF49VX48kubfBYmAkkEm4G6fq9jfdvyqeoWVe2uqk2Bwb5tuwFEpCrwITBYVRcW9gGqmqGqyaqaXLu23TkyJuT17u2GmO7b50YUffih1xGZUxBIIlgE1BeReiJSAegFTPc/QERqiUjeuR4EJvq2VwDexXUkvxW8sI0xnktJcZ3I9evDNdfAiBHWiVxGFZkIVPUQMACYBawG3lDVlSIyXESu9R3WDlgjIt8BZwHpvu2/B9oAfURkqe/RJNjfhDHGI7Gx8NlnbtLZAw+4TuSDB72Oypwkq1lsjDl1ubmuxsHjj8MVV7i6yNWrex1VRLGaxcYYb+V1Ir/4omshXHYZrFvndVQmQJYIjDHBc8st8MknsH07NG8OX3zhdUQmAJYIjDHB1aYNLFzobg1dcQW89prXEZkiWCIwxgRf/fpuJnLz5nDjjTB8uI0oCmGWCIwxJaNmTfj4Y3e7aOhQ6NvXlqUIUeW8DsAYE8ZOOw0mT3ZrFQ0bBps3uxFFZ5zhdWTGj7UIjDElS8S1CCZNgnnzoHVrpozdSUKCG2yUkABTpngcY4SzFoExpnT06QOxsUzp8ir976lEjq/LYMMG6N/fPU9L8yy6iGYtAmNM6enYkcE1xpGjMUdtzsmBwYM9islYIjDGlK6NP1UofPvGUg7E5LNEYIwpVXFxx9le14aXesUSgTGmVKWnQ8zRd4aI4VfSz58Ihw8X/iZToiwRGGNKVVoaZGRAfLwbUBQfp2Rc8z5pn/Zzq5ju2+d1iBHHEoExptSlpcH69W7R0vUbhLTpvWDsWJg+HTp2hF27vA4xolgiMMaEhgED4I03IDMTWrVy40pNqbBEYIwJHddf75al2LrVLWW9fLnXEUWEgBKBiKSKyBoRWSsigwrZHy8ic0RkuYjME5FYv323isj3vsetwQzeGBOG2rZ1NQ3ArWSa99yUmCITgYhEA88CnYBEoLeIJBY4bASuLnESMBx4wvfeGsBQoDmQAgwVEStbZIw5sYsucquXnn02XHUVvP++1xGFtUBaBCnAWlVdp6oHgKlA1wLHJAKf+p7P9dv/f8BsVd2lqj8f9mveAAAPlElEQVQDs4HUUw/bGBP24uLgv/91SaFbN7d4nSkRgSSCOsAmv9dZvm3+lgHdfc+7AVVEpGaA7zXGmMLVqgVz5kD79m4Z6xEjvI4oLAWrs/h+oK2ILAHaApuBgGeGiEh/EckUkczs7OwghWSMCQtVqsAHH8Dvfw8PPAB/+5sVuQmyQFYf3QzU9Xsd69uWT1W34GsRiEhloIeq7haRzUC7Au+dV/ADVDUDyABITk62f2FjzNFOOw1efdUVu3nqKcjOdrPSytkCysEQSItgEVBfROqJSAWgFzDd/wARqSUieed6EJjoez4LuEpEqvs6ia/ybTPGmJMTHQ3PPnuktsH118P+/V5HFRaKTASqeggYgPsFvhp4Q1VXishwEbnWd1g7YI2IfAecBaT73rsLeAyXTBYBw33bjDHm5Im4Smdjx8J778HVV8PevV5HVeaJhti9tuTkZM3MzPQ6DGNMqHv5ZdeBnJwMM2ZAjRpeR+QpEVmsqsnFea/NLDbGlE033+zqHy9Z4iahbd3qdURlliUCY0zZdd11rjXw44/QurVbyc6cNEsExpiyrUMH+OQTt2Lp5ZfD6tVeR1TmWCIwxpR9LVrAf/7jCtu0bg2LF3sdUZliicAYEx4uusgtUFe5MlxxBXz+udcRlRmWCIwx4eP8810yOOsst1jdp58W/R5jicAYE2bq1oX58+G886BzZ9eZbE7IEoExJvycfTbMmweNGrmRRW+/7XVEIc0SgTEmPNWs6VYuTU6Gnj1hyhSvIwpZlgiMMeGrWjVX+rJNGzcBbcIEryMKSZYIjDHhrXJl+PBDSE2F/v3dOkXmKJYIjDHhr1IlmDbN9RfcfTeMHOl1RCHFEoExJjJUqABvvAE33AD33QdPPul1RCHDqjoYYyJH+fKuwE358vDgg3DgADzyiNdRec4SgTEmspQrBy+95L4OHQoHD8Lw4a7WQYSyRGCMiTzR0a7KWYUK8PjjrmXw5JMRmwwC6iMQkVQRWSMia0VkUCH740RkrogsEZHlItLZt728iLwoIitEZLWIPBjsb8AYY4olKgrGj4c77nB1kO+7D0KsUFdpKbJFICLRwLPAlUAWsEhEpqvqKr/DhuBKWD4vIonADCABuAE4TVUvEpEYYJWIvKaq64P8fRhjzMmLinJ1kCtUgFGj3Oqlo0dHXMsgkFtDKcBaVV0HICJTga6AfyJQoKrv+RnAFr/tp4tIOaAScAD4XxDiNsaY4BBxSSAqyn2FiEsGgSSCOsAmv9dZQPMCxwwDPhaRu4DTgY6+7W/hksZWIAa414rXG2NCjgj885/ua94cgwhKBsGaR9AbmKyqsUBn4GURicK1Jg4D5wL1gPtE5LyCbxaR/iKSKSKZ2dnZQQrJGGNOggiMGAEDB8KYMUxJfYmEeCUqChISwnupokBaBJuBun6vY33b/N0GpAKo6gIRqQjUAm4EZqrqQWC7iHwOJAPr/N+sqhlABkBycnJk9tYYY7znSwZTvr2E/jOuIwfXItiwwa1OAZCW5mF8JSSQFsEioL6I1BORCkAvYHqBYzYCHQBEpCFQEcj2bb/Ct/10oAXwbXBCN8aYEiDC4G9uJIfTj9qckwODB3sUUwkrMhGo6iFgADALWI0bHbRSRIaLyLW+w+4D/igiy4DXgD6qqrjRRpVFZCUuoUxS1eUl8Y0YY0ywbNxUeN/Axo2lHEgpCWhCmarOwA0J9d/2iN/zVUCrQt63FzeE1Bhjyoy4OHc76JjtdRUIvw5kW3TOGGMKSE+HmJijt8XwK+mNpoTlpDNLBMYYU0BaGmRkQHy86z+Oj1MyrnqLtI9udovVhVkysLWGjDGmEGlp/iOEBPQWuPNL+Mc/3IJ1jz0WNvMMLBEYY0wgROBf/4JDh9y9o/Ll3eqlYcASgTHGBCoqCsaNc8lg2DDXMgiDMaWWCIwx5mRERcGECS4ZDBniksHf/uZ1VKfEEoExxpysvHoGhw7BoEFQsSLcc4/XURWbJQJjjCmO6GhX6ey33+Avf4FKlY6sQ1HG2PBRY4wprnLl4LXXoHNn+NOfXGIogywRGGPMqahQAd5+G664Avr2hddf9zqik2aJwBhjTlXFivDee9CqlZt88N57Xkd0UiwRGGNMMJx+OnzwASQnw+9/DzNneh1RwCwRGGNMsFStCh99BImJ0K0bzJ3rdUQBsURgjDHBVL06zJ4Nv/sdXHMNLFzodURFskRgjDHBVquWSwZnnw2dOsGyZV5HdEKWCIwxpiSccw7MmQOVK8NVV8GaNV5HdFwBJQIRSRWRNSKyVkQGFbI/TkTmisgSEVkuIp399iWJyAIRWSkiK3z1jI0xJvzFx7tkANCxY+HVbkJAkYlARKJxJSc7AYlAbxFJLHDYEFwJy6a4msbP+d5bDngF+JOqNgLaAQeDFr0xxoS6Bg3g449h717o0AG2bvU6omME0iJIAdaq6jpVPQBMBboWOEaBqr7nZwBbfM+vApar6jIAVd2pqodPPWxjjClDLr7YjSb66Se48krYudPriI4SSCKoA2zye53l2+ZvGHCTiGThahvf5dveAFARmSUiX4vIX08xXmOMKZtatID334e1ayE1Ff73P68jyheszuLewGRVjQU6Ay+LSBRuUbvLgTTf124i0qHgm0Wkv4hkikhmdnZ2kEIyxpgQ0769W45i6VK49lrYt8/riIDAEsFmoK7f61jfNn+3AW8AqOoCoCJQC9d6mK+qO1Q1B9dauKTgB6hqhqomq2py7dq1T/67MMaYsuLqq+Hll2H+fOjZEw56320aSCJYBNQXkXoiUgHXGTy9wDEbgQ4AItIQlwiygVnARSIS4+s4bgusClbwxhhTJvXqBc89524V/eEPkJvraThF1iNQ1UMiMgD3Sz0amKiqK0VkOJCpqtOB+4AJInIvruO4j6oq8LOIjMQlEwVmqOqHJfXNGGNMmfGnP8HPP8NDD7nZyM884+oieyCgwjSqOgN3W8d/2yN+z1cBrY7z3ldwQ0iNMcb4GzQIdu2CESOgRg1XB9kDVqHMGGO8IgJPPeVaBo8+6loGHpS8tERgjDFeEoHx42H3blfysnp1uOWWUg3BEoExxngtOhqmTIE9e1zn8RlnQNeC83ZLji06Z4wxoeC00+Ddd6FZMzeqaP78UvtoSwTGGBMqKleGDz+EhAQ34ayUlq+2RGCMMaGkVi2YNQuqVHFLUaxbV+IfaYnAGGNCTVycSwYHDrhaBtu2lejHWSIwxphQlJjobhNt3eqqnJXgInWWCIwxJlS1aOEWqVuxAq67DvbvL5GPsURgjDGhLDUVXnwR5s6FtDQ4HPySLpYIjDEm1N14I4weDe+8A3feCapBPb1NKDPGmLLgnntchbMnn4Q6deDhh4N2aksExhhTVvz977BlCzzyCJxzDvTrF5TTWiIwxpiyQgT+/W/Yvh1uvx3OPNNNPDtF1kdgjDFlSfny8OabbimKnj3hiy9O+ZSWCIwxpqzJW4oiNpYpV04i4dwDuMxQPAElAhFJFZE1IrJWRAYVsj9OROaKyBIRWS4inQvZv1dE7i9uoMYYY/zUrs2UP/+X/jnPsGFrhVM6VZGJQESigWeBTkAi0FtEEgscNgR4Q1Wb4moaP1dg/0jgo1OK1BhjzFEGP3MWOcSc8nkCaRGkAGtVdZ2qHgCmAgUXylagqu/5GcCWvB0ich3wI7DylKM1xhiTb+PG4JwnkERQB9jk9zrLt83fMOAmEcnC1Ta+C0BEKgN/Ax495UiNMcYcJS4uOOcJVmdxb2CyqsYCnYGXRSQKlyBGqereE71ZRPqLSKaIZGZnZwcpJGOMCW/p6RBz6neGAppHsBmo6/c61rfN321AKoCqLhCRikAtoDlwvYg8BVQDckVkv6r+y//NqpoBZAAkJycHd+60McaEqbQ093XwYNiwofjnCaRFsAioLyL1RKQCrjN4eoFjNgIdAESkIVARyFbV1qqaoKoJwGjg7wWTgDHGmOJLS4P16wEWLy7uOYpMBKp6CBgAzAJW40YHrRSR4SKSN6XtPuCPIrIMeA3ooxrkVZGMMcaUCAm139fJycmamZnpdRjGGFOmiMhiVU0uznttZrExxkQ4SwTGGBPhLBEYY0yEC7k+AhH5BVjjdRwhohaww+sgQoRdiyPsWhxh1+KIC1S1SnHeGIr1CNYUt8Mj3IhIpl0Lx67FEXYtjrBrcYSIFHuUjd0aMsaYCGeJwBhjIlwoJoIMrwMIIXYtjrBrcYRdiyPsWhxR7GsRcp3FxhhjSlcotgiMMcaUIs8SQQDlL08Tkdd9+78UkYTSj7J0BHAtBorIKl8Z0DkiEu9FnKWhqGvhd1wPEVERCdsRI4FcCxH5ve9nY6WIvFraMZaWUy2XGy5EZKKIbBeRb46zX0RkjO86LReRSwI6saqW+gOIBn4AzgMqAMuAxALH/BkY53veC3jdi1hD5Fq0B2J8z++I5GvhO64KMB9YCCR7HbeHPxf1gSVAdd/rM72O28NrkQHc4XueCKz3Ou4SuhZtgEuAb46zvzOuLLAALYAvAzmvVy2CQMpfdgVe9D1/C+ggIlKKMZaWIq+Fqs5V1Rzfy4W4mhDhKJCfC4DHgH8A+0szuFIWyLX4I/Csqv4MoKrbSznG0nJK5XLDiarOB3ad4JCuwEvqLASqicg5RZ3Xq0QQSPnL/GPULYW9B6hZKtGVrkCuhb/bcBk/HBV5LXxN3bqq+mFpBuaBQH4uGgANRORzEVkoIqmlFl3pKna53Ah0sr9PgNCcWWyOQ0RuApKBtl7H4gVf+dORQB+PQwkV5XC3h9rhWonzReQiVd3taVTeyCuX+08RaYkrl9tYVXO9Dqws8KpFEEj5y/xjRKQcrrm3s1SiK12BXAtEpCMwGLhWVX8rpdhKW1HXogrQGJgnIutx90Cnh2mHcSA/F1nAdFU9qKo/At/hEkO4CbRc7hvgyuXiqiTWKpXoQktAv08K8ioRBFL+cjpwq+/59cCn6usNCTNFXgsRaQqMxyWBcL0PDEVcC1Xdo6q19Ej504W4axKOlYwC+T8yDdcaQERq4W4VrSvNIEtJscvllmqUoWE6cItv9FALYI+qbi3qTZ7cGlLVQyKSV/4yGpiovvKXQKaqTgdewDXv1uI6R3p5EWtJC/BaPA1UBt709ZdvVNVrj3vSMirAaxERArwWs4CrRGQVcBh4QFXDrtUc4LW4D5ggIvfiOo77hOMfjiLyGi751/L1hwwFygOo6jhc/0hnYC2QA/QN6LxheK2MMcacBJtZbIwxEc4SgTHGRDhLBMYYE+EsERhjTISzRGCMMRHOEoExxkQ4SwTGGBPhLBEYY0yE+3/ByMCUOYIixwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "def f(x):\n",
    "    if not isinstance(x, float):\n",
    "        x = float(x)\n",
    "    return math.sin(x) / x\n",
    "\n",
    "# actual result (by sympy): Si(1) = 0.946083070367183\n",
    "\n",
    "I, m, p = adaptsim(f, 1e-64, 1)\n",
    "\n",
    "# 画图\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "xf = np.linspace(1e-64, 1)\n",
    "yf = [f(x) for x in xf]\n",
    "\n",
    "yp = [f(x) for x in p]\n",
    "\n",
    "plt.plot(xf, yf, 'r-', label='f=sin(x) / x')\n",
    "plt.plot(p, yp, 'bo', label='adaptive points')\n",
    "plt.legend()\n",
    "plt.xlim(0, 1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "I=-1.4260192969841916\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4VOX1wPHvm7AGENmJBBKQVSCAYVWxIEuxqLhV9BdbQGxcULQoqMQqWiPuWrQWgyJY0traulURChUsLiARN0BEShbCGkBACIGQnN8f7yQkzIRMMpO5s5zP88wzM/feufdwuZkz7303IyIopZRSUU4HoJRSKjhoQlBKKQVoQlBKKeWiCUEppRSgCUEppZSLJgSllFKAJgSllFIumhCUUkoBmhCUUkq51HE6gPJatmwpCQkJToehlFIh5YsvvtgrIq183U9QJYSEhAQyMzOdDkMppUKKMSbHH/vRW0ZKKaUATQhKKaVcNCEopZQCgqwOQalQUFRURF5eHoWFhU6HoiJMgwYNiIuLo27durWyf00ISlVTXl4eTZo0ISEhAWOM0+GoCCEi7Nu3j7y8PDp27Fgrx9BbRkpVU2FhIS1atNBkoALKGEOLFi1qtWSqCaEaMjIgIQGiouxzRobTESmnaDJQTqjt605vGXkpY2ERKTcZCo7ZU5aTAykpdl1ysoOBKaWUn2gJ4XSKi2H5cpg8mdRJO8uSQamCAkhNdSg2FdHmzJlDjx49SPby18gNN9xA69at6dWrV4Xl+/fvZ9SoUXTp0oVRo0bx448/AvZ+9dSpU+ncuTOJiYmsW7eu0n3ffPPNfPLJJ17HvmDBArKzsyk/n/upy6ZPn0737t1JTEzkiiuu4MCBA277Kd121qxZFd5X5cYbb2Tjxo1Vbvfcc8/x2muvAfDGG2/Qs2dPoqKi3DrPzp49m86dO9OtWzeWLl1a6f4ee+wxMry8rfDCCy8wf/58r7b1KxEJmkdSUpIEhTVrRKZOFWnbVgREmjQRQ4mAuD2McTpYFWgbN250OgTp1q2bbNu2zevtP/roI/niiy+kZ8+eFZZPnz5dZs+eLSIis2fPlhkzZoiIyPvvvy9jxoyRkpIS+eyzz2TgwIGV7rtPnz5y4sSJKmPIy8uTyZMny8MPPyx//vOfJSUlxeMyEZGlS5dKUVGRiIjMmDGjLK7ynnnmGXn55Zfl7rvvlpkzZ8rSpUu9OxleKCoqkt69e5fFsHHjRtm0aZP87Gc/k7Vr15Ztt2HDBklMTJTCwkLZunWrdOrUqdJzMWzYMNmzZ49Xxz9y5Ij07dvX4zpP1x+QKX74DvbLFzlwJvAPYBPwHTAEaA4sA35wPTeraj9BkRBeesmelvr1Ra68UuSNN0QKCiQ+3j0ZgEh8y8NOR6wCzOmEcNNNN0ndunWlV69e8swzz3j9uaysLLeE0LVrV9mxY4eIiOzYsUO6du0qIiIpKSnyl7/8xeN25W3cuFF++ctfui2/7LLLZOHChSIiMnfuXPm///s/ERHZtWuXxMfHy9ixY6W4uLjSZeW9+eabZZ8/1ezZs6VevXry3//+123d4cOH5Re/+IUkJiZKz5495fXXXxcRqfCl3qhRI5k5c6YkJibKoEGDZNeuXSJiE9KECRPc9nlqQnj00Ufl0UcfLXs/evRo+fTTT90+d/DgQTnvvPPclk+dOlUeeughERFZsmSJDB06tOwcXH755bJmzRq3z9RmQvBXHcIfgCUicrUxph4QA8wE/iMijxlj7gXuBe7x0/H8KiPD3vrJzRU6mItJ6/YQyWvugKZNy7ZJS7N1BgUFJz8XE1VI2t7fwNSW8NRTUK+eA9ErR915J3z1lX/32bcvPPdcpavnzp3LkiVLWLFiBS1btmTFihX89re/ddsuJiaGTz/99LSH2r17N7GxsQC0bduW3bt3A7B9+3bat29ftl1cXBzbt28v27bUBx98wJgxY9z2m56ezvnnn0/Hjh15+umnWb16NTt27ODBBx/khhtuoGPHjkyZMoXf/e53bsv+9Kc/VdjX/PnzGT9+vNsx/vCHP9CqVSumTp3KkiVLKCwsZNSoUWXrlyxZwllnncX7778PwMGDB932ceTIEQYPHkxaWhozZsxg3rx53H///XzyySckJSWd9tyVnqfBgwe7nadTLV++nBEjRrgtnz17NgMGDGDo0KFMnTqVxYsXExVl7+T379+fVatWMXDgwCrj8Bef6xCMMU2BC4FXAETkuIgcAMYBC12bLQQu9/VYtSEjw37R5+SAiCGnpD0pOalkvNe0wnbJyZCeDvHxYIx9Tn+1DsnTYuH552HYMPBwIShV24YPH85XX33l9qgqGZzKGFPtVixLly71mBDatGnDww8/zPDhw3n66adp3rw5Z511FvPmzaNDhw4MHTqUF1980eOy8tLS0qhTp47HupKpU6cyefJkGjVqRFpaGiNHjqywvnfv3ixbtox77rmHVatW0bRpU7d91KtXj0suuQSApKQksrOzAdi5cyetWvk8eGiZJUuWcPHFF7stj4mJYd68eYwaNYrbbruNs88+u2xd69at2bFjh99i8IY/SggdgXzgVWNMH+AL4A6gjYjsdG2zC2jjh2P5XWpqxV/9AAWF0aSmurceSk4+dVkd4GkYNAhuuAFGjoQvv4QGDWo5ahU0TvNLPlB8KSG0adOGnTt3Ehsby86dO2ndujUA7dq1Y9u2bWXb5eXl0a5duwqfLSgo4MCBA5x11lke9/3tt9/SokULty+1iRMnum3radmCBQt47733+M9//uMxUZUuK61UPnWbrl27sm7dOhYvXsz999/PiBEjeOCBBypsU7du3bLPRUdHc+LECQAaNmzoVXt/b84TwOeff+5W8ilV2XkqLCykYcOGVcbgT/5oZVQHOBf4k4j0A45gbw+Vcd3j8tgEwBiTYozJNMZk5ufn+yGc6snNrd5yj665Bt56CzZtgt//3i9xKeUtX0oIl112GQsX2oL8woULGTduXNny1157DRFh9erVNG3a1O120YoVKxg+fLjH/X7++ed88MEHfPnllzz11FNkZWVV69+0ZMkSnnjiCd59911iYmKq9dlSO3bsICYmhuuvv57p06eftqXUqXr06MGWLVuq3O6yyy7j9ddf59ixY2RlZfHDDz+43eLZsGED3bt3Jzo62u3zOTk5PP3003z55Zd88MEHrFmzpmzd5s2b3VqF1TpfKyGAtkB2ufdDgfeB74FY17JY4Puq9uVEpXKllcXxNdjZpEki0dEiX37p5yhVMHG6UllEJD4+XvLz873e/tprr5W2bdtKnTp1pF27dvLyyy+LiMjevXvloosuks6dO8uIESNk3759IiJSUlIit956q3Tq1El69epVoSK11JQpU2TFihVuywsLCyUxMVG++OILERF55513ZNiwYVJSUuJ1vGeffbbExcVJnz59pE+fPnLTTTd5/dlSS5Yskd69e0ufPn2kf//+Zf+GUyuVS73xxhtlFcnZ2dkydOjQsnVvvvmmtGvXTurVqyetW7eW0aNHl6175JFHpFOnTtK1a1dZvHixWxxPPvmkvPrqq27LS0pKZMSIEfLOO++IiEhmZqb06tVLjh49KiIi/fr1k71797p9LhRaGa0CurlezwKedD3udS27F3iiqv04kRAWLRKJiS6skAxiYuzyatu/X6RNG5F+/URczdVU+AmGhBAM+vXrJ8ePH3c6jFpz+eWXy+bNm33ez8iRIz220DqddevWyfXXX+9xXW0mBH91TLsdyDDGfAP0BR4FHgNGGWN+AEa63ged5MRvSS++gfgzfjxZWZxew97HzZrBH/8IX35JxvWLdZgLFdbWrVtXa6NuBoPHHnuMnTt3Vr1hFZYtW+Z2u60qe/fu5fcO3H42NrkEh/79+0vAp9C89FJYtQq2boXmzf2yy4wBz5KSmUIBjcqWxcT4kGhUUPnuu+/o0aOH02GoCOXp+jPGfCEi/X3dd2QPXfHxx/Dee3DPPX5LBgCpO2+vkAxAh7lQSgW/yE0IInDvvRAbC3fc4ddd5+7w3Jq3Wi2XlFIqwCI3Ibz3HnzyCTzwgL2f40cdOlRvuVJKBYPITAjFxXDffdClC0ye7Pfdp6W555iYGCEtze+HUkopv4nMhPDhh7BhA8yaBbXQSqLCMBcI8WSTPvEzrVBWAbVgwQJuu+22Gn02Ozubv/zlL2XvMzMzmTp1qr9Cq5bzzjuvym2ee+45Ck4dckBVW0QlhLIZz0aPJIFsMgqvrLVjJSdDdjaUnCghu9+VJH9wPRQV1drxVPAKxZn2Tk0I/fv3Z86cOY7E4k2Pa00I/hExCaHCIHYYcogn5fYGtf/HGR0NDz8MWVnw5z/X8sFUsKk4eOLJmfZ8ve4uv/xykpKS6NmzJ+np6WXLX331Vbp27crAgQMrTFrzr3/9i0GDBtGvXz9GjhxZNqrprFmz+NWvfsWQIUPo0qUL8+bNA+Dee+9l1apV9O3bl2effZaVK1dyySWXUFJSQkJCQoUJa7p06cLu3bvJz8/nqquuYsCAAQwYMMDjpDkLFixg3LhxDBs2jC5duvDQQw+VrXvmmWfo1asXvXr14rlyY0Q1btwYgJUrVzJs2DCuvvpqunfvTnJyMiLCnDlz2LFjB8OHD2f48OEUFxczceJEevXqRe/evXn22Wd9O9mRxB+92/z1qM2eyn4doqK6SkpEkpJEOnYUCeOenZGiOj2Va+u6Kx1ioqCgQHr27Cl79+6VHTt2SPv27WXPnj1y7NgxOe+882TKlCkiIrJ///6yoSPmzZsn06ZNExGRBx98UBITE6WgoEDy8/MlLi5Otm/fLitWrJCxY8eWHa/8+6lTp8r8+fNFRGT16tUyYsQIERG57rrrZNWqVSIikpOTI927d3eL+9VXX5W2bdvK3r17y2Jfu3Zt2bANhw8flp9++knOOeccWbdunYicHF5ixYoVcsYZZ8i2bdukuLhYBg8eXHa88kN5ZGZmysiRI8uO+eOPP/p2soNMKMyHEPT8MohdTRlj6ysuvRQWLYJJkwJwUBUMauu6mzNnDm+99RYA27Zt44cffmDXrl0MGzasbNjm8ePHs3nzZsCOwjl+/Hh27tzJ8ePH6dixY9m+xo0bR8OGDWnYsCHDhw/n888/58wzz6z02OPHj+fhhx9m0qRJvP7662VzFSxfvrzC1JSHDh3i8OHDZb/wS40aNYoWLVoAcOWVV/Lxxx9jjOGKK66gUaNGZctXrVpFv379Knx24MCBxMXFAdC3b1+ys7O54IILKmzTqVMntm7dyu23387YsWMZPXq0l2dVRcwtI8ebgo4dC0lJ8MgjWpcQQWrjulu5ciXLly/ns88+4+uvv6Zfv35VDtV8++23c9ttt/Htt9/y0ksvVdj+1GGjq5oTYciQIWzZsoX8/HzefvttrrzS1sWVlJSwevXqstFWt2/f7pYManK88urXr1/2uvxw1eU1a9aMr7/+mmHDhjF37lxuvPFGr/cf6SImIXhuCkrgmoIaAw8+aIfIWLQoQAdVTquN6+7gwYM0a9aMmJgYNm3axOrVqwEYNGgQH330Efv27aOoqIg33nijwmdKx+kvHe661DvvvENhYSH79u1j5cqVDBgwgCZNmvDTTz95PH7pr/lp06bRo0ePsl/7o0eP5vnnny/b7qtKZpJbtmwZ+/fv5+jRo7z99tucf/75DB06lLfffpuCggKOHDnCW2+9xdChQ70+J+Xj3bt3LyUlJVx11VU88sgj1Rr2OtJFzC2j5GRgz25Spx0ll3g6xBvS0gI8ttAll8C559pSwvXX10qTVxVcSq8vO0WrLRn4et2NGTOGuXPn0qNHD7p161Y2hWNsbCyzZs1iyJAhnHnmmfTt27fsM7NmzeKXv/wlzZo146KLLqowP0FiYiLDhw9n7969/O53v+Oss86iVatWREdH06dPHyZOnOh262b8+PEMGDCABQsWlC2bM2cOU6ZMITExkRMnTnDhhRcyd+5ct/gHDhzIVVddRV5eHtdffz39+9sheCZOnFg2l8CNN97odszTSUlJYcyYMZx11lk899xzTJo0iZKSEsBOU6m8E1mD2730Etx8s53Iplu32jvO6fzrX3DZZTB/vtYlhKhwGtxu1qxZNG7cmLvvvjsgx1uwYAGZmZm88MILATleONLB7fxl+XKIi4OuXZ2LoXwpQesSlFJBJHJKCMXF0Lq1/XX+6qu1cwxvlZYSXn0VPMwlq4JbOJUQVOjREoI/fPUV7N8PI0c6HUnFUoKHVhIq+AXTDykVOWr7uouchLB8uX0OhoRQ2uLof//TFkchqEGDBuzbt0+TggooEWHfvn00aNCg1o4RObeMRo2C3bvhm29qZ//VJWL7JRw6ZCu560RMg6+QV1RURF5eXpVt/5XytwYNGhAXF+c2dam/bhlFxrfQ0aN2msxbb3U6kpNKey+PG2cHtpkwwemIlJfq1q1boaevUuEiMm4ZffIJHDsWHLeLyrv0UjLiZpDwm1FERUnIjISplApPkVFCWL7cdgK78EKnI6kg4y+GlD1pFBTZ/4bSkTAhwB3mlFKKSCkhLF8OQ4aAh3FVnJSaCgXHK+bkggK7XCmlAi38E8K+fbBuXfDdLsLhEViVUuoU4Z8QVqywLXqCMCE4PgKrUkqVE/4JYflyaNIEBgxwOhI3HkfCbFgSuBFYlVKqHL8lBGNMtDHmS2PMe673HY0xa4wxW4wxfzPG1PPXsbxRNo/tSy+ScGILGX8Lvvrz5GRIT4f4eDBGiCeH9JF/1wplpZQj/FlCuAP4rtz7x4FnRaQz8CMw2Y/HOq2K8ydHkXO0tV/msa0NycmQnQ0lJYbsq+4iedUtcPiw02EppSKQXxKCMSYOGAu87HpvgIuAf7g2WQhc7o9jeSM11bbWKS8kWu/cdRccOOD84HtKqYjkrxLCc8AMoMT1vgVwQERKR27LA9r56VhVCtnWO0OGwHnnwbPP2tFZlVIqgHxOCMaYS4A9IvJFDT+fYozJNMZk5ufn+xoOEOKtd+66C7KywDWBulJKBYo/SgjnA5cZY7KB17G3iv4AnGmMKa3JjQO2e/qwiKSLSH8R6d+qVSs/hONqvdOgpMKygM6f7Itx4+Dss+Gpp2xzWaWUChCfE4KI3CcicSKSAFwLfCgiycAK4GrXZhOAd3w9lreSkyH9lnXEk21b78Tb1jwh0XonOhp++1tYswY+/dTpaJRSEaQ2+yHcA0wzxmzB1im8UovHcpPcYinZdKTk4GGys0MkGZSaOBGaNYOnn3Y6EqVUBPFr43wRWQmsdL3eCgz05/6rZf162xGhSRPHQqixRo3glltg9mzYsgU6d3Y6IqVUBAjfnsobNkCvXk5HUXO33WYnzfnDH5yORCkVIcIzIRQV2VnIevZ0OpKai42F666zfRIOHHA6GqVUBAjPhPDDDzYphHIJAWzl8pEjMG+e05EopSJAeCaEDRvsc6gnhL59YfhwmDPHJjillKpF4ZkQ1q+HqCjo3t3pSHz3299CXh78859OR6KUCnPhmxA6d4YGDZyOxHdjx0KXLvDMM9pRTSlVq8IzIYR6C6PyoqLgzjth7VrtqKaUqlXhlxAKC22lcii3MDrVhAm2o9qzzzodiVIqjIVfQti0CUpKwqeEANCoERkXvEjCP58mKkpISAjOuR2UUqEt+KYR81W4tDAqJyMDUpZfQwFRIHbin5QUuy6khuRQSgW18CshrF8PdevaitgwkZoKBUcr/leFxIQ/SqmQEp4JoVs3mxTCRMhO+KOUCinhlxDCqYWRS0hP+KOUChnhlRAOH7azjYVTCyNcE/7EVFwWU+9EaEz4o5QKGeGVEDZutM9hVkJITrYT/MTHYyf8qbOd9PhHtEJZKeVX4ZUQwrCFUankZMjOhpISQ/bT/yT5h4fgs8+cDkspFUbCKyGsX2+Hq+jY0elIatcNN9iOak8+6XQkSqkwEn4J4Zxz7LzE4axxY5gyBd5+GzZvdjoapVSYCK+EEIYtjCp1++1Qr57Ou6yU8pvwSQg//gjbt4ddC6NKtW4NEyfCwoWwa5fT0SilwkBYJISMDEjo0ZAoikl46rbIGefnrrvg+HF4/nmnI1FKhYGQTwgZGXZcn5zdDRCiyMmPISUlQgZ/69IFrrwSXnwRfvrJ6WiUUiEu5BNCaqod16e8iBrnZ/p0OHAAXnnF6UiUUiEu5BNCxI/zM2gQXHihnVFN511WSvkg5BOCjvMDzJgB27bB3/7mdCRKqRDmc0IwxrQ3xqwwxmw0xmwwxtzhWt7cGLPMGPOD67mZ7+G6S0uDmIYV5xqOiSGyxvm5+GLbuuqJJ3TeZaVUjfmjhHACuEtEzgEGA1OMMecA9wL/EZEuwH9c7/0uORnS79xIPNl2nJ94O+5PRI3zExVl6xK+/Rb+/W+no1FKhSifE4KI7BSRda7XPwHfAe2AccBC12YLgct9PVZlkmM/JJuOlOzYTXZ2hCWDUtddB+3a2VKCUkrVgF/rEIwxCUA/YA3QRkR2ulbtAtpU8pkUY0ymMSYzPz+/Zgdevx5atIA2Hg8RGerVgzvvhA8/hMxMp6NRSoUgvyUEY0xj4J/AnSJyqPw6ERHA481tEUkXkf4i0r9Vq1Y1O/j69fYeujE1+3y4SEmBM87QQe+UUjXil4RgjKmLTQYZIvKma/FuY0ysa30ssMcfx3IjElljGJ3OGWfArbfCP/4BmzY5HY1SKsT4o5WRAV4BvhORZ8qteheY4Ho9AXjH12N5tH07HDyoCaHUtGlkRP+ahKQWREVBQkKE9NpWSvmsjh/2cT7wK+BbY8xXrmUzgceAvxtjJgM5wDV+OJa70klxImVQuypk/LsVKfISBQX1AMjJsXeSIEIr25VSXvM5IYjIx0BlN+9H+Lr/Kq1fb581IQCuoTxO1KuwrHQoD00ISqnTCfmeyqxfD7GxtpWR0qE8lFI1FvoJYcMGLR2Uo0N5KKVqKrQTggh8/z107+50JEEjLc0O3VFeTIOSyBrKQylVI6GdEA4cgEOHoFMnpyMJGsnJduiO+HjsUB4ml/T+L2n9gVKqSqGdELKy7HNCgqNhBJvkZMjOhpISQ/Zv/0DyZ7fDli1Oh6WUCnKhnRCys+1zx46OhhHUpk+HunXhkUecjkQpFeRCOyFoCaFqbdvCLbfAokVaSlBKnVZoJ4TsbDjzTPtQlZsxw5YStGZZKXUaoZ0QsrK0dOCNtm3h5pvhz3+G//3P6WiUUkEq9BOC1h94p7SU8PvfOx2JUipIhW5CELG3jLSE4J3YWDsS6muvwTffOB2NUioIhW5CyM+3g/RoCcF7qam2vuXuu3XuZaWUm9BNCKVNTrWE4L3mzeGBB2DZMli61OlolFJBJnQTQmmTUy0hVM+tt8LZZ9tSwokTTkejlAoioZsQtIRQM/XqweOP20EB5893OhqlVBAJ3YSQlQUtW0Ljxk5HEnquvBIuuMDePvrpJ6ejUUoFidBNCNrCqOaMgaeegt274YknnI5GKRUkQjchaB8E3wwaRMbgOSSk/YaoKNG5l5VSIZoQSkrsZMFaQqixjAxI+XoKOdIBEVM297ImBaUiV2gmhF274NgxLSH4IDUVCo5W/O8vnXtZKRWZQjMh6LDXPtO5l5VSpwrNhKDDXvtM515WSp0qtBNCfLyzcYQwj3MvmwLSHih0JiCllONCLiFkZEDCYzcTRTEJPRpqJWgNVZx7GeLbFJIuN5K88X6nQ1NKOaTWE4IxZowx5ntjzBZjzL2+7Csjw7aEyTnSEiFKW8b46OTcy5C9qwHJNzeFZ5+FtWudDk0p5QAjtTjqpTEmGtgMjALygLXAdSKy0dP2/fv3l8zMzEr3l5BgW5ueKj7+ZD2z8sHBg3DOObYHeGamnT9BKRX0jDFfiEh/X/dT2yWEgcAWEdkqIseB14FxNd2ZtoypZU2bwosv2vkSnnrK6WiUUgFW2wmhHbCt3Ps817Ia0ZYxATBuHFx9NTz0EGze7HQ0SqkAcrxS2RiTYozJNMZk5ufnn3bbtDSIqV9cYVlMjM4d73fPPw8NG9oKmpISp6NRSlUiI6O09X1Skj/2V9sJYTvQvtz7ONeyMiKSLiL9RaR/q1atTruz5GRIn/AJ8WRjjBAfb1vKJCf7P/CI1rYtPP00fPQRvPyy09EopTwoa2TjoV61pmo7IawFuhhjOhpj6gHXAu/6ssPk2A/JNp0oKSwiO1uTQa2ZNAkuugimT4cdO5yORilV3rFjpM4ooqDAv7ut49/dVSQiJ4wxtwFLgWhgvohs8Gmn2dnQrp2d6EXVHmPgpZegd2+47TZ4802nI1IqshQV2U64mzfbxw8/wJYt9pGbS25Jkd8PWasJAUBEFgOL/bZDHfY6cDp3tpXL99wDf/sbjB/vdERKhZ/9++G77+xj0yb4/nubALZurTjNbfPm0KULnH8+TJhAhxcKyNnn3wnCaj0h+F12NvzsZ05HETmmTYM33yRjwlJSp11J7s66dOhgK/L1dp1SXhKBPXvs1LXr159MAN99Z5eXatDAfun37g1XXQXdukHXrvbRokWFXaZ1sXUI/rxtFFoJoagI8vK0hBBIdeqQ8X/vkbImhoIdtqNaaQ9x0KSglJtDh2xfnvXrTyaA9eth796T25x5JvToAZdcYp9LH/HxEB3t1WFK//ZSU/1XsVyrPZWrq6qeymzdCmefbSeHnzQpcIFFOO0hriJNRob9os3NpfIScUmJ/U76+mv7+OYb+1z+j6JxY+jVyz569jz53LatrafzE3/1VA6tEkLpKKdaQggo7SGuIklpc87SWzG2RCyQk0Ny7ApYtw6+/NJ++R8+bDeKirK3dQYNgt/8BhIT7W2fDh38+sVf20IzIeg8CAHVoYPnEoL2EFfhKHWmUFBQ8Uu8oMCQmgrJ3ACNGkHfvjBxIvTrZ7/8e/a0nTlDXGglhOxse38tLs7pSCJKWpp75VWMOUraAwZo4FhcSvmsuNhW7K5dawd0XLuW3NzVgPuv+lwTD5u+t63vohwf5KFWhFZCyMqyP0vrhFbYoa585VVuLnRodZS0PTeSvDIaJi0MqSKxinC7dsGaNScfa9fCTz/ZdU2aQFISHc44SM6hZm4f7dAf1aNNAAAShklEQVTB2NtCYSy0vlmzs/V2kUOSk8tXqjWEh7rCrFm2CfDkyQ5GplQliopsRe+nn558lFZ81aljb/X86lf2vv/AgfbLPiqKtAwPJeIIGTMttBJCVhZcfLHTUSiA+++Hjz+2vZgTE2HAAKcjUpHuxx8rfvl//vnJb/W4OBgyBO64wyaAc8+t9J6/W4k4gvrdhE5CKCyEnTu1hBAsoqNtc4zBg22S/vhj6N7d6ahUJMnLg1Wr7OPjj21bfxH7679fP9va57zzbCJo377q/ZVTsUQcOUInIZQ2c9Emp8GjdWv497/hggtg9Gj45JNq/+Ep5RUR2+Z/5Uo7Cu+qVSfb+zdubL/4r7nGXosDB9p7PKraQichlP7nawkhuHTuDEuW2LqE0aPtH2rLlk5HpUKdiB3ErTQBrFwJ210j57dqBRdeaG//DB0KffpoQxM/CZ2zqJ3SglffvvDuu/Dzn8PYsfCf/9hfbUqd4rQ9gHNy4MMPTz5Kh11v0waGDbM/OoYNs7cmtWVbrQidhJCdbYe8jo11OhLlyc9+Bn//O1x5pX38619Qv77TUakg4rEH8A0n4OVXSc59zN4SAnsr8qKL7Jf/sGG29Y8mgIAInd4VWVl28Jww7RASFi67zM6wtmwZ/PrXttOPUi6p95W4jcxZcLwOqR/93LZUmzPHVgzv2gV//SvcdJMd7VOTQcCEVglB6w+C38SJdlTH6dPt+O0vvqh/0JHqxAnb+3fZMli+nNxtKzxulkt7eOutAAenPAmdhJCVBVdc4XQUyht3322TwuOP216gr7yit48iRU6ObXm2dKmtSzpwwP4gOPfc0/cAVkEhNBLC4cOQn68lhFAye7YdCuD++23rkDffhGbuXwYqxB0+bFsALV1qE8HmzXZ5XJytSxo9GkaMgJYtI7oHcKgIjRvy2gch9Bhjm5MsWmT7J5x/vk6eEGIyMuxvsKgo+5yRgW0O+s038OST9ou+RQu49FI7R0nnzvDcc7Bxo21G9MordtpVVzPk5GRIT7dVgcbY5/T0yOwAFqxCY4Kc996zF91nn9mesSq0rFxpb/fVrw/vvw9JSU5HpKpwaosggJjoY6Q3mUbygRftgt69bVPjMWNshzC9LegYf02QExolhNJfllpCCE3DhtlSQoMGtkPRe+85HZE6nZISUu8+5t4iqLg+qcdn2V/+eXkVSwqaDMJCaCSErCw7EFXr1k5HomrqnHNg9Wo7b+y4cfDCC2RkiPstCeWM/Hz7H3D99dCmDbm76nrcLPdoK7jhBmjXLsABqkAIjUrl0ian2nwxtLVta28fXXcdGbd/Skr0jRQU2wl27DSFdjO9pxwAxcV2PoAlS+wjM9PWD7RsCT//OR2WFJCzz723uc6SF95CIyFs26ZXYrho3BjeeYfUlocp+LHibGsFBbYeWhNCLdmxw7YGWrLE9g348UdbPBs8GB56yNYFJCVF/JwAkSw0EkJenu3JqMJDVBS5B87wuKp0/hJVPR7HCLr6mB0WujQJfPut3bhtW3vb7uKLYeRI24HwFJE8J0Ak8ykhGGOeBC4FjgP/AyaJyAHXuvuAyUAxMFVEltboIEVFtiu7zqMcVjp0ONmauMLyFkeARgGPJ5R5HCNoQiFMuoXkogVQt64dFfTxx20poHdvr26/RuqcAJHM10rlZUAvEUkENgP3ARhjzgGuBXoCY4AXjTHRNTrCzp323qZWYoWVtDT3IetjzFHS9v7Gdmhav96ZwEJNfj6pdx7x0CKoAan1n7KDDO7fb3sNz5hhS9paF6cq4VNCEJF/i8gJ19vVQOnP+HHA6yJyTESygC3AwBodpHQMdC0hhBWPnZQW1CU5rRcsX25/xV5zjSaGUx06ZPty3HWXnRWsdWty93qeCjL3SAu45BIdilx5zZ/NTm8APnC9bgdsK7cuz7XMjTEmxRiTaYzJzM/Pd98gL88+a0IIO8nJtgFZSYl9Tv51HZg50zYzTk2FDz6wv2jHj4cNG5wO1xkFBfbXfWqqnQqyeXP7Jf/HP9rXjzxCh7ZFHj+q7TBUdVWZEIwxy40x6z08xpXbJhU4AVS7JbmIpItIfxHp36pVK/cNNCFEnhYt4JFHbJaYORMWL7YlhvHjbVPJIOpdX1Meh4UAe3vn3XftaLGDB0PTprbi9/HH7cb33WcnjzlwoCxRpD1V3/32m7YIUjUhIj49gInAZ0BMuWX3AfeVe78UGFLVvpKSksTNtGkiDRuKlJS4r1ORYe9ekZkzRRo3FgGRTp3s+/Xr3TZdtEgkPl7EGPu8aFHAo63SokUiMTH2n1L6iKlTKIvaTT+5oF49kQsuELnvPpHFi0UOHapyn8H+71a1B8gUH7/LxV55PiWDMcBGoNUpy3sCXwP1gY7AViC6qv15TAjjx4t06eKn06ZC2oEDIvPni4waJRIVZS/f3r1FHn1UZOtWz1+0MUHy5VhSIpKVJfLGGxJ/xv4KMZY+4hvsEklLE/nvf0WOHnU6YhVC/JUQfBrczhizxfWlv8+1aLWI3Oxal4qtVzgB3CkiH3jey0keB7e74AI7deaHH9Y4ThWGdu+2U3b+9a920EMgoU4eOSfcq6ri430baPW08wCfSgT27LGV4eUfGzbYuSGAKIoRD3drjbH1KUpVl78Gtwv+0U4TEuyAaK+95khMKgRkZ8M//kHU9Gmev2gRSua8YJsux8XZ57ZtIbrqltAeR/1sWEL6/bkkn/OVPXZOzsnnrCx7f79Uixa2/qNXL/tISiLhqiRyct2bfvqauFTkCvuEkJEBqTOF3FyhwxkHSXuxmXaSUaeVkOC5s1s8OWSTUHFhVBTExkKbNrbjVp06NkGUfy4uJmHlAnKOx3rYZzbZuEbfbdTIHjw+3j66dTuZAFq3dmv37zHJxOjcAKrm/JUQgnLoipN/MAYw5BxqpgOfqSqlpVUy/s5L7WH0btunpfSRl2ef9+yxveGLi+0cwMeP2+fiYjvExvE2Ho+Va+JhbaZNBM2bV6uzlw4LoYJVUJYQKv2lp0VqVYVq3e/3gl6LKhSE9QQ5lQ1wpgOfqaq4dXbz8Ve3xyE2tI2/ClNBmRAq62GpPS9VoOk8wCqSBGVC0F9lKpj4u9ShVLAKyoRQ9qus0V4MJfqrTCmlAiAoWxmBayz2V66xrT4+/tjpcJRSKuwFZQmhTF6ezoOglFIBErwJQcS2E9dRTpVSKiCCNyEcOGB7GGlCUEqpgAjehKDzICilVEAFf0LQOgSllAqI4E0IOpeyUkoFVPAmhLw82zU01n2kSaWUUv4X3AmhbVs7NLFSSqlaF9wJQesPlFIqYII3IWgfBKWUCqjgTQh5eZoQlFIqgIIzIRw+bDumaUJQSqmACc6EoE1OlVIq4II7IWilslJKBUxwJgQdtkIppQIuuBOClhCUUipg/JIQjDF3GWPEGNPS9d4YY+YYY7YYY74xxpxbrR3m5UGLFtCwoT/CU0op5QWfE4Ixpj0wGsgtt/hioIvrkQL8qVo73b5dSwdKKRVg/ighPAvMAKTcsnHAa2KtBs40xng/KJH2QVBKqYDzKSEYY8YB20Xk61NWtQO2lXuf51rmHU0ISikVcHWq2sAYsxxo62FVKjATe7uoxowxKdjbSnTo0AGOHYM9ezQhKKVUgFWZEERkpKflxpjeQEfga2MMQBywzhgzENgOtC+3eZxrmaf9pwPpAP379xd27rQrtA5BKaUCqsa3jETkWxFpLSIJIpKAvS10rojsAt4Ffu1qbTQYOCgiO73asfZBUEopR1RZQqihxcAvgC1AATDJ609qQlBKKUf4LSG4SgmlrwWYUqMdaUJQSilHBF9P5e3boXFjOOMMpyNRSqmIEnwJQZucKqWUIzQhKKWUAjQhKKWUcgm+hLBzp/ZBUEopBwRXQigqguJiLSEopZQDgishHD9unzUhKKVUwAVXQigqss+aEJRSKuCCKyGUlhC0DkEppQIu+BJCvXrQsqXTkSilVMQJroRQVGRvF9nRU5VSSgVQcCWE48e1/kAppRwSfAlB6w+UUsoRwZUQSm8ZKaWUCrjgSggimhCUUsohwZUQQBOCUko5JKgSwhckkTBlLBkZTkeilFKRJ6gSAkDOrvqkpKBJQSmlAizoEgJAQQGkpjodhVJKRZagTAgAublOR6CUUpElaBNChw5OR6CUUpElKBNCTAykpTkdhVJKRZagSwjx8ZCeDsnJTkeilFKRpY7TAZSXlASZmU5HoZRSkSnoSghKKaWc4XNCMMbcbozZZIzZYIx5otzy+4wxW4wx3xtjfu7rcZRSStUun24ZGWOGA+OAPiJyzBjT2rX8HOBaoCdwFrDcGNNVRIp9DVgppVTt8LWEcAvwmIgcAxCRPa7l44DXReSYiGQBW4CBPh5LKaVULfI1IXQFhhpj1hhjPjLGDHAtbwdsK7ddnmuZG2NMijEm0xiTmZ+f72M4SimlaqrKW0bGmOVAWw+rUl2fbw4MBgYAfzfGdKpOACKSDqS7jvWTMeb76nzeIS2BvU4H4QWN079CIc5QiBE0Tn/r5o+dVJkQRGRkZeuMMbcAb4qIAJ8bY0qwJ3A70L7cpnGuZVX5XkT6e7Gdo4wxmRqn/2ic/hMKMYLG6W/GGL802Pf1ltHbwHAAY0xXoB42m74LXGuMqW+M6Qh0AT738VhKKaVqka8d0+YD840x64HjwARXaWGDMebvwEbgBDBFWxgppVRw8ykhiMhx4PpK1qUB1R2RKN2XeAJI4/QvjdN/QiFG0Dj9zS9xGvuDXimlVKTToSuUUkoBAUoIxpj5xpg9rroGT+uNMWaOa6iLb4wx55ZbN8EY84PrMcHhOJNd8X1rjPnUGNOn3Lps1/Kv/FXj70Ocw4wxB12xfGWMeaDcujGu4US2GGPudTjO6eViXG+MKTbGNHetC8j5NMa0N8asMMZsdA2/coeHbRy/Pr2M0/Hr08s4Hb8+vYwzGK7PBsaYz40xX7vifMjDNvWNMX9znbM1xpiEcuuqN4SQiNT6A7gQOBdYX8n6XwAfAAbbp2GNa3lzYKvruZnrdTMH4zyv9PjAxaVxut5nAy2D5HwOA97zsDwa+B/QCdsi7GvgHKfiPGXbS4EPA30+gVjgXNfrJsDmU89JMFyfXsbp+PXpZZyOX5/exBkk16cBGrte1wXWAINP2eZWYK7r9bXA31yvz3Gdw/pAR9e5jT7d8QJSQhCR/wL7T7PJOOA1sVYDZxpjYoGfA8tEZL+I/AgsA8Y4FaeIfOqKA2A1tn9FwHlxPiszENgiIlvFNgh4HXvua0U147wO+GttxVIZEdkpIutcr38CvsO9V73j16c3cQbD9enl+axMwK7PGsTp1PUpInLY9bau63Fqxe84YKHr9T+AEcYYQw2GEAqWOoTKhrrweggMB0zG/mosJcC/jTFfGGNSHIqpvCGuYuYHxpiermVBeT6NMTHYL9J/llsc8PPpKmr3w/4KKy+ors/TxFme49dnFXEGzfVZ1fl0+vo0xkQbY74C9mB/gFR6fYrICeAg0IIanM+gmiAnVBg7yutk4IJyiy8Qke3Gjvi6zBizyfUL2QnrgHgROWyM+QW2A2EXh2LxxqXAJyJSvjQR0PNpjGmM/YO/U0QO1dZxfOVNnMFwfVYRZ9Bcn17+vzt6fYrtw9XXGHMm8JYxppeIeKyX81WwlBAqG+qipkNg1BpjTCLwMjBORPaVLheR7a7nPcBbODi6q4gcKi1mishioK4xxpchRWrbtZxSHA/k+TTG1MV+KWSIyJseNgmK69OLOIPi+qwqzmC5Pr05ny6OXp/ljnkAWIH7bcmy82aMqQM0BfZRk/NZ25Ui5So+Eqi8EnQsFSvtPnctbw5kYSvsmrleN3cwzg7Y+3DnnbK8EdCk3OtPgTEOxtmWk31MBgK5rnNbB1vx2ZGTlXY9nYrTtb4ptp6hkRPn03VeXgOeO802jl+fXsbp+PXpZZyOX5/exBkk12cr4EzX64bAKuCSU7aZQsVK5b+7XvekYqXyVqqoVA7ILSNjzF+xLQtaGmPygAexlSOIyFxgMbYlxxagAJjkWrffGPN7YK1rVw9LxWJboON8AHtv7kVbZ8MJsQNftcEW5cBe1H8RkSUOxnk1cIsx5gRwFLhW7BVywhhzG7AU26JjvohscDBOgCuAf4vIkXIfDeT5PB/4FfCt6z4twEzsl2swXZ/exBkM16c3cQbD9elNnOD89RkLLDTGRGPv6PxdRN4zxjwMZIrIu8ArwJ+NMVuwyeta17+h2kMIaU9lpZRSQPDUISillHKYJgSllFKAJgSllFIumhCUUkoBmhCUUkq5aEJQSikFaEJQSinloglBKaUUAP8PdzZzB/o+1GkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def f(x):\n",
    "    return (100 / x**2) * math.sin(10 / x)\n",
    "\n",
    "I, m, p = adaptsim(f, 1, 3, eps=1e-4)\n",
    "print(f'I={I}')\n",
    "\n",
    "xf = np.linspace(1, 3)\n",
    "yf = [f(x) for x in xf]\n",
    "\n",
    "yp = [f(x) for x in p]\n",
    "\n",
    "plt.plot(xf, yf, 'r-', label='f=100 / x**2 * sin(10 / x)')\n",
    "plt.plot(p, yp, 'bo', label='adaptive points')\n",
    "plt.legend()\n",
    "plt.xlim(1, 3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt4FOXZP/DvnZCQLIlgEiAEyCYROQWCSASqRFEERZGIVUGCxUOlFaj0IFiNCsUGqFZe4EfFFxUBsxat0mqRwwtIi6KAAUE5CkKCOYA5ABJIQpK9f3/M7rKH2WSTPc3u3p/r2ovdmdmZZ5fN3DPP4X6ImSGEEEKE+bsAQgghtEECghBCCAASEIQQQphIQBBCCAFAAoIQQggTCQhCCCEASEAQQghhIgFBCCEEAAkIQgghTNr4uwDWEhISOCUlxd/FEEKIgLJnz54KZu7o7n40FRBSUlJQUFDg72IIIURAIaIiT+xHqoyEEEIAkIAghBDCRAKCEEIIABprQxAiENTX16O4uBi1tbX+LooIMVFRUejWrRsiIiK8sn8JCEK0UHFxMWJjY5GSkgIi8ndxRIhgZlRWVqK4uBipqaleOYZUGTXDYABSUoCwMOVfg8HfJRL+Vltbi/j4eAkGwqeICPHx8V69M5WA4ITBACQkAJMmAUVFALPy76RJAJHySEiQABGqJBgIf/D2704CggqDAZgyBaisbHq7ykolQEhgEEIEAwkIVszVQ5MmAZcuuf6+ykrg4YeBqVO9VjQhbCxZsgR9+vRBTk6OS9s/9thj6NSpE/r162ezvKqqCiNHjsS1116LkSNH4uzZswCU+uqnnnoKPXr0QEZGBvbu3et037/+9a+xY8cOl8u+cuVKFBYWwno+d/tlM2fORO/evZGRkYFx48bh3LlzDvsxbztnzhyb18355S9/iUOHDjW73aJFi7B69WoAwD/+8Q+kp6cjLCzMYfDs/Pnz0aNHD/Tq1QubNm1yur8FCxbA4OKV49KlS7FixQqXtvUoZtbMY9CgQewv+fnMOh2zUjnUugeRsh8R3A4dOuTvInCvXr34hx9+cHn7//73v7xnzx5OT0+3WT5z5kyeP38+MzPPnz+fZ82axczMn3zyCd95551sNBr5yy+/5MGDBzvd94ABA7ihoaHZMhQXF/Pjjz/Oc+fO5XfeeYenTJmiuoyZedOmTVxfX8/MzLNmzbKUy9rChQv5zTff5Keffpqfe+453rRpk2tfhgvq6+u5f//+ljIcOnSIjxw5wrfccgt/9dVXlu0OHjzIGRkZXFtbyydOnOC0tDSn38Xw4cP5xx9/dOn4Fy9e5Ouuu051ndrvD0ABe+Ac7JETOYAOAD4AcATAYQA/AxAHYDOAY6Z/r25uP/4MCHq9e8HA/IiP99tHED7i74Dwq1/9iiMiIrhfv368cOFCl9938uRJh4DQs2dPLi0tZWbm0tJS7tmzJzMzT5kyhd99913V7awdOnSIH3jgAYflY8eO5VWrVjEz8+uvv84TJ05kZubTp0+zXq/nu+++mxsbG50us7Z27VrL++3Nnz+fIyMjefv27Q7rqqur+a677uKMjAxOT0/nNWvWMDPbnNTbtWvHzz33HGdkZPCQIUP49OnTzKwEpMmTJzvs0z4gzJs3j+fNm2d5PWrUKP7iiy8c3nf+/Hm+8cYbHZY/9dRT/Kc//YmZmTdu3MhZWVmW7+Dee+/lXbt2ObzHmwHBU91OFwPYyMz3E1EkAB2A5wBsZeYFRPRHAH8E8IyHjudRBoPSYOwJlZVKm8LixYCLd/MikP32t8C+fZ7d53XXAYsWOV39+uuvY+PGjdi2bRsSEhKwbds2/O53v3PYTqfT4YsvvmjyUGfOnEGXLl0AAImJiThz5gwAoKSkBN27d7ds161bN5SUlFi2NduwYQPuvPNOh/0uX74cN910E1JTU/Hqq69i586dKC0txezZs/HYY48hNTUV06ZNwwsvvOCwbNmyZTb7WrFiBcaPH+9wjMWLF6Njx4546qmnsHHjRtTW1mLkyJGW9Rs3bkRSUhI++eQTAMD58+cd9nHx4kUMHToUeXl5mDVrFt544w08//zz2LFjBwYNGtTkd2f+noYOHerwPdnbsmULRowY4bB8/vz5uOGGG5CVlYWnnnoK69evR1iYUpOfmZmJzz77DIMHD262HJ7idhsCEbUHcDOAtwCAmS8z8zkA2QBWmTZbBeBed4/lDeYG5Kbo9UB+vvKIj29+n5WVyj6loVn4wq233op9+/Y5PJoLBvaIqMW9WDZt2qQaEDp37oy5c+fi1ltvxauvvoq4uDgkJSXhjTfeQHJyMrKysvDaa6+pLrOWl5eHNm3aqLaVPPXUU3j88cfRrl075OXl4fbbb7dZ379/f2zevBnPPPMMPvvsM7Rv395hH5GRkRgzZgwAYNCgQSgsLAQAlJWVoWNHt5OHWmzcuBGjR492WK7T6fDGG29g5MiRmD59Oq655hrLuk6dOqG0tNRjZXCFJ+4QUgGUA3ibiAYA2ANgBoDOzFxm2uY0gM4eOJbH5eY6b0DW6YDly22v9M3Pp04F7C5kbFy6BMyYIXcJQa+JK3lfcecOoXPnzigrK0OXLl1QVlaGTp06AQC6du2KH374wbJdcXExunbtavPeS5cu4dy5c0hKSlLd97fffov4+HiHk9ojjzzisK3aspUrV2LdunXYunWraqAyLzM3Kttv07NnT+zduxfr16/H888/jxEjRuDFF1+02SYiIsLyvvDwcDQ0NAAAoqOjXerv78r3BAC7d+92uPMxc/Y91dbWIjo6utkyeJInehm1AXA9gGXMPBDARSjVQxamOi7VLgBENIWICoiooLy83APFaZlTp5yvsw8G1l57rfm7hcpKuUsQ3ufOHcLYsWOxapVyI79q1SpkZ2dblq9evRrMjJ07d6J9+/YO1UXbtm3Drbfeqrrf3bt3Y8OGDfj666/x17/+FSdPnmzRZ9q4cSNefvllfPzxx9DpdC16r1lpaSl0Oh0mTZqEmTNnNtlTyl6fPn1w/PjxZrcbO3Ys1qxZg7q6Opw8eRLHjh1zqOI5ePAgevfujfDwcIf3FxUV4dVXX8XXX3+NDRs2YNeuXZZ13333nUOvMK9ztxECQCKAQqvXWQA+AXAUQBfTsi4Ajja3L380KjtrTNbrm3+vKz2TXNmPCCz+blRmZtbr9VxeXu7y9hMmTODExERu06YNd+3ald98801mZq6oqODbbruNe/TowSNGjODKykpmZjYajTx16lROS0vjfv362TSkmk2bNo23bdvmsLy2tpYzMjJ4z549zMz80Ucf8fDhw9loNLpc3muuuYa7devGAwYM4AEDBvCvfvUrl99rtnHjRu7fvz8PGDCAMzMzLZ/BvlHZ7B//+IelIbmwsJCzsrIs69auXctdu3blyMhI7tSpE48aNcqy7s9//jOnpaVxz549ef369Q7leOWVV/jtt992WG40GnnEiBH80UcfMTNzQUEB9+vXj2tqapiZeeDAgVxRUeHwvkDoZfQZgF6m53MAvGJ6/NG07I8AXm5uP/4ICPn5zLrwWpuTuE7nevfR/HylZ1FTXVFFcNFCQNCCgQMH8uXLl/1dDK+59957+bvvvnN7P7fffrtqD62m7N27lydNmqS6zpsBwVMD034DwEBE3wC4DsA8AAsAjCSiYwBuN73WDPMgtIcfZkQ3/oT4qIsgUhqQm6oqspeTA1RUOK8+YpYcSCI47d2712tZN7VgwYIFKCsra37DZmzevNmhuq05FRUVeOmll9w+dkt5pNspM+8DkKmyyrGflQaYexYpjcmESnSEjhjvvNP6RuDFi633aauo6EpPJmlkFiIw9OrVC7169fLLsa27z/pSSKauUOtZdKmGkJvb+n3m5Ch3Fno9oNZ+bu51JIKDcpcuhG95+3cXkgHB2SC0pnocuSInBygsdJ6RUHodBYeoqChUVlZKUBA+xazMhxAVFeW1Y4TcBDkGg5K6Wu1vOTnZM8dITnYedHJzpdoo0HXr1g3FxcXwRzdpEdrMM6Z5S8gFhNxc9WBABOTleeYYeXlKxlQ1nkqRIfwnIiLCazNWCeFPIVdl5KxaiNlzV+45Oc57HRFJtZEQQptCKiAYDMpUmGqUxmDPWbxYOfnbY4ZbjddCCOEtIRMQzF1NGxsd1+l0nqsuMsvJUa+aAtxvvBZCCG8ImYDgLIldeHjLBqK1hLO7Dk81Xgth7cCBA7jjjjuQkJDgkbl39+3bh0GDBkGn02HQoEHYZ5Xme/To0YiJibE8IiMj0b9/f7ePKfwrZAKCs6tyo9F7vX7y8pS7D2sERlGRjF4WnhcREYEHH3wQb731ltv7unz5MrKzszFp0iScPXsWkydPRnZ2Ni5fvgxAmQehurra8rjxxhvxwAMPuH1c4WeeyH/hqYc3cxk5yzfk7eRz+fnmBHpGJjS2OmeSCBwlJSV83333cUJCAqekpPDixYuZmXn27Nn885//nB988EGOiYnhgQMH8r59+yzvW7BgASclJXFMTAz37NmTt2zZ0qrjHzt2jJU/bdfKpWbTpk2clJRkk5Cue/fuvGHDBodtT548yWFhYXzy5MlWlVe4DxrLZaRpBgPw00+OyyMjPd92YM88WE2vJ9h/3ZcuSQNzsDEajbjnnnswYMAAlJSUYOvWrVi0aJFl8vWPPvoIDzzwAKqqqjBx4kTce++9qK+vx9GjR7F06VJ89dVXuHDhAjZt2oSUlBQAwLvvvosOHTo4fZxyoVGquXLZO3jwIDIyMmyqnjIyMnDw4EGHbVevXo2srCxLeUXgComAkJsL1Nc7Lo+N9d0gMWd/s9LAHFy++uorlJeX48UXX0RkZCTS0tLwxBNPYM2aNQCUWbnuv/9+RERE4Pe//z1qa2uxc+dOhIeHo66uDocOHUJ9fT1SUlIss2dNnDgR586dc/pIdqFRqrly2auurnaYYax9+/a4cOGCw7arV69WneBGBJ6QGJjm7KRbVeW7MjgbvSwNzMGlqKgIpaWl6NChg2VZY2MjsrKyoNfrbeYpDgsLQ7du3VBaWoqsrCwsWrQIc+bMwcGDB3HHHXdg4cKFTmcj82S5ACAmJsay/NChQ4iJicFPdrfVP/30E2JjY22Wff755zh9+jTuv/9+j5RT+FdI3CE4O+n68mSs1sAMMKqrpXE5mHTv3h2pqak2V/AXLlzA+vXrAcBmukWj0Yji4mLLSX/ixIn4/PPPUVRUBCLCM888AwAwGAw2PXrsH65UGTVXLusG4uTkZKSnp+Obb76xydf0zTffID093Wa/q1atwn333WcTUETgComAcNddAJHtoABvjD1oijkbqjKC2VwWQmWlMj5CgkJwGDx4MGJjY/GXv/wFNTU1aGxsxIEDB/DVV18BAPbs2YO1a9eioaEBixYtQtu2bTF06FAcPXoUn376Kerq6hAVFYXo6GiEmUZR5uTk2Jyw7R/mKiNmRm1traUnUG1tLerq6lwql73hw4cjPDwcS5YsQV1dHZYuXQoAuO222yzb1NTU4P3335fqomDiiZZpTz280ctIbZpLIuYnn/T4oVzizpSdIjCUlJTwhAkTuHPnztyhQwceMmQIb9682aGX0XXXXWeZZnL//v18ww03cExMDF999dV89913c0lJSYuOe/LkSfPc5ZaH3uqH5axczuzdu5evv/56joqK4oEDB/LevXtt1r/77rucnJzcoqkxhXfAQ72MiJ0Np/WDzMxMLigo8Og+U1LU6+71eqX3j6+FhTlPrmc0+r48wnfmzJmD48ePIz8/399FEUGGiPYws9okZS0S9FVGWuvdo4X2DCGEUBP0AUFrJ2C1xmVdZINP2zOEEEJNUAcEgwGorgbsp7T0dYOyNeupNokY8WFnEd14AQ8/zJLOIsjNmTNHqouEpgVtQDBnN62sBADzaEtGfLz3ktm5yjx6+Z13CDVtYlHZeDWYCUVF0uNICOE/QRsQ1LObEmJitDOFZW4ucOmy7djAS5eAGTP8VCDhNzExMThx4oS/iyFCXNAGBK01JqtxVpbKSrlLCDXV1dVIS0tzaVsiwvHjxz16/KVLlyIzMxNt27ZVHVewdetW9O7dGzqdDrfeeiuKVLruVVVVoWPHjhg2bJhHyyZ8J2gDgtYak9U0VRZJeid8KSkpCc8//zwee+wxh3UVFRW477778NJLL6GqqgqZmZkYP368w3bPPPMM+vTp44viCi8J2oCg2pvHj43JapoqS1GR3CUEopSUFMyfPx99+/bF1VdfjUcffRS1tbUAgDfeeAM9evRAXFwcxo4di9LSUsv7rK/6H3nkEUybNg133303YmNjMWTIEHz//fcAgJtvvhkAMGDAAMTExOC9995DRUUFxowZgw4dOiAuLg5ZWVkwtnBQy3333Yd7770X8SqTga9duxbp6el44IEHEBUVhTlz5mD//v04cuSIZZsvvvgCBw4cwKOPPtqyL0xoiscCAhGFE9HXRLTO9DqViHYR0XEieo+IIj11LFdFRwPmQZtaaEy2l5NjTmWhThqYA5PBYMCmTZvw/fff47vvvsOf//xnfPrpp3j22Wfx/vvvo6ysDHq9HhMmTHC6jzVr1mD27Nk4e/YsevTogVzTLeP27dsBAPv370d1dTXGjx+PV199Fd26dUN5eTnOnDmDefPmWdJWmwOF2mPMmDEufZ6DBw9iwIABltft2rXDNddcY0mF3djYiOnTp2Pp0qUemalN+I8n7xBmADhs9fovAP6HmXsAOAvgcQ8eq0mOPYwINTW+OnrLLF6slvROIfMlBKbp06eje/fuiIuLQ25uLv7+97/DYDDgsccew/XXX4+2bdti/vz5+PLLL1HoZLj8uHHjMHjwYLRp0wY5OTk201fai4iIQFlZGYqKihAREYGsrCzLiXndunVO02avW7fOpc/TXCrsJUuWYMiQIRg0aJBL+xPa5ZGAQETdANwN4E3TawJwG4APTJusAnCvJ47lCrUeRlo9uZrHJTijpUZw4RrrFNd6vR6lpaUoLS2F3mqS7ZiYGMTHx6OkpER1H4mJiZbnOp0O1cqAGlUzZ85Ejx49MGrUKKSlpWHBggUe+BRXNJUKu7S0FEuWLEGelupiRat56g5hEYBZAMwVl/EAzjFzg+l1MYCuHjpWswKhh5G1nBxloJoaLTWCC9dYp7g+deoUkpKSkJSUZNMz5+LFi6isrETXru7/WcTGxuLVV1/FiRMn8PHHH2PhwoXYunUrAGD06NFO02aPHj3apf2np6dj//79NmX//vvvkZ6ejt27d6OsrAx9+/ZFYmIiZsyYgd27dyMxMRGNjY1ufzbhW24HBCIaA+BHZt7TyvdPIaICIiooLy93tzgAAqOHkb1AaAQXrvnb3/6G4uJiVFVVIS8vD+PHj8dDDz2Et99+G/v27UNdXR2ee+45DBkypFXTTnbu3NlmzMK6detw/PhxMDPat2+P8PBwS+rsDRs2OE2bvWHDBss+GhoaUFtbi8bGRjQ2NqK2thYNDcr13Lhx43DgwAF8+OGHqK2txdy5c5GRkYHevXtj9OjRKCwsxL59+7Bv3z7MnTsXAwcOxL59+xAeHu7eFyl8z910qQDmQ7kDKARwGsAlAAYAFQDamLb5GYBNze3LU+mv8/OZdVGBN6F9fj6zPtnIQCOHo54BI+v12i+3uEKv1/O8efO4T58+3L59e/7FL37BFy9eZGbmZcuWcVpamiW99Q8//GB5HwA+duwYMzNPnjyZc3NzLeu2bdvGXbt2tbxetmwZJyYmcvv27fm9997jhQsXsl6vZ51Ox127duW5c+e2uNyzZ892SJ09e/Zsy/rNmzdzr169OCoqim+55RY+efKk6n7efvttvummm1p8fOEeeCj9tUfnMwAwHMA60/N/AJhgev46gKnNvd+TASE+tpYBIwNGjo8PnJNqfj6zLuJywAUzodDr9U3OMSCEN3gqIHhzHMIzAH5PRMehtCm85cVjWVh6GF1oC633MFKTmwtcqo+wWabVBnEhRHDxaEBg5v8w8xjT8xPMPJiZezDzA8xc58ljORNIPYzUBFqDuBAieLRpfpPAEugn1ORk9RnetNwgLq5wNq5AiEAQdKkrArGHkTXV3kbRLL2NhBBeF3QBIS8P0EXZ5nEJpO6bDhPooBzRVIOHH4ZMoCOE8KqgCwg5OcDyJ3ZDj0IQMfR67eUwao7NBDphsai8pAMzZAKdAJCSkoItW7a4tY8XXngB/fv3R5s2bTBnzhyX3nP58mX06dMH3bp1syz77rvvkJ2djY4dOyIuLg533HEHjh49avO+EydOYMyYMYiNjUVCQgJmzZplWVdVVYVx48ahXbt20Ov1ePfddy3rPvnkEwwbNgwdOnRAYmIifvnLX1pSWVhTS4ltMBhsBsjpdDoQEfbsUYYyzZkzBxERETbbmMddNPeZVq5cifDwcJv3/uc//3HpOxRBGBAMBiB3dW+cQjKSuylVLYEUDKzl5gKXjFE2ywKpgVy0To8ePfDyyy/j7rvvdvk9r7zyCjp27Giz7Ny5cxg7diyOHj2KM2fOYPDgwcjOzrasv3z5MkaOHInbbrsNp0+fRnFxMSZNmmRZP23aNERGRuLMmTMwGAx48sknLQntzp8/j+effx6lpaU4fPgwSkpKMHPmTIdyqaXEzsnJsRkg99prryEtLQ3XX3+9ZZvx48fbbGOeK6K5zwQAP/vZz2zeO3z4cJe/x5Dnib6rnnq4Ow4hP1/psx8sffiJbD+L+UHk75IJZzw5DiEnJ8dmcJgzJ06c4N69e/P69ettBrDZq6ysZABcUVHBzMz/+7//y8OGDVPdtrq6miMiIvjo0aOWZZMmTeJnnnlGdfsPP/yQ+/XrZ7Nsx44dPHToUF6xYkWTg9WGDx/Oc+bMsbyePXs25+TkON2+qc8UqgPjEADjEHwu0Luc2gv0BvJQ9vnnnztNO92hQwd8/vnnHjvWb37zG8ybNw/RSr53p7Zv347ExETLnAc7d+5ESkoKRo8ejYSEBAwfPhzffvstAKVqpk2bNujZs6fl/QMGDLDcIajtOz093fLa1ZTYRUVF2L59O37xi1/YLP/3v/+NuLg4pKenY9myZS5/JgD4+uuvkZCQgJ49e+Kll16ypOAQzQuqbqeB3uXUXl6e0mZgHeQCqYE8lA0bNgznzp3z+nH++c9/orGxEePGjWuyrry4uBjTpk3DwoULbZZt27YNH3/8MUaMGIHFixcjOzsbR44cQXV1Na666iqbfVinvLa2efNmrFq1Crt27bIss06JbQ4yalavXo2srCykpqZalj344IOYMmUKOnfujF27duHnP/85OnTogIceeqjZz3TzzTfjwIED0Ov1OHjwIMaPH482bdrg2WefdVoGcUVQ3SEE2xW1dY8jgBGOBly6xMjNlYblYJGenm5p/Pzss89a9N6LFy9i1qxZWLJkSZPblZeXY9SoUZg6darNSTU6OhrDhg3D6NGjERkZiaeffhqVlZU4fPhwkymvre3cuRMTJ07EBx98YLmbaElK7NWrV2Py5Mk2y/r27YukpCSEh4fjxhtvxIwZM/DBBx/YbOPsM6WlpSE1NRVhYWHo378/XnzxRYf3CueC6g4hLw+Y8lgDLl2+8rEC/Yra3CA+5QnGpRrlc5l7G1mvF9ry2WefNZleesOGDcjKynJaBeOKY8eOobCwEFlZWQCURuLz588jMTHRUh109uxZjBo1CmPHjrXMumaWkZGBHTt2qO67Z8+eaGhowLFjx3DttdcCUGZps64W+vrrrzF27FisWLECI0aMsCy3TokNADU1NaipqUFiYiJKSkosWVB37NiB0tJS3H///U1+TiIy50oDgCY/U3PvFc3wREOEpx6eSG6Xf8/fWY9CJgqeTKF6vXrjsl7v75IJe55oVL58+TLX1NTwQw89xLm5uVxTU8MNDQ0O29XX13NZWZnl8eGHH3KXLl24rKyMGxoa+Pz583zDDTfwtGnTVI9z5MgRjo6O5s2bN3NDQwMvXLiQ09LSuK6ujpmZx48fzxMmTODq6mr+/PPP+aqrruIDBw4wM/O3337LnTp14jVr1jjst7a21qZcixYt4sGDB3NZWZnNdk888QQ//PDDDu//17/+xVVVVWw0GnnXrl2clJTEK1euZGZu9jOtX7+eT58+zczMhw8f5vT0dJsG62AFLWY7dffhkWynY8Yw2/V2CHTS2yhweCIgTJ482SEV9dtvv83MzNu3b+d27dqpvs8+TfbKlSsZAOt0Om7Xrp3lUVRUZNnmww8/5GuuuYZjY2P5lltusZzwmZUePNnZ2azT6bh79+5sMBgs6x555BEmIpv99u3bV7Vcaj1/ampquH379rxlyxaH7SdMmMBxcXHcrl077tWrFy9evNjlz/SHP/yBO3XqxDqdjlNTU/mFF17gy5cvO/uqg4anAgKxhm6nMjMzuaCgoFXvNRiU3kSnioxI1lUib3nHoKlOSUlRz2+k1ysD2IQQoY2I9jBzprv7CYpGZXPK66IigBGGoksdg2pEr2p+o8iGgG4bEUJoT1AEhGAbf2DPIb8RVSHaWC35jYQQHhUUASHYxh+osclv1CYWlQ0dwJLfSAjhQUEREIJt/EFTZEY1IYS3BEVAUK1jD/DxB86Ewt2QEMI/giIgWOrYw4tBCMyU164KpbshIYRvBUVAAICcu86isLE7jH95BYWFwRkMAGczqhmD8m5ICOFbQREQDAYgpU80wtCIlL9OD+oGVtUZ1SAzqgkh3BfwAcEyBuFMlDIGoVwX9L1ubHochceisqad9DgSQrgt4Ecqh/Io3lD+7EKIK2Skskko97oJ5c8uhPC8gA8IodzrJpQ/uxDC89wOCETUnYi2EdEhIjpIRDNMy+OIaDMRHTP9e7X7xXWUlwfoom2rvYJ1DII9tR5HAKO6WtoRhBAt54k7hAYAf2DmvgCGAphGRH0B/BHAVma+FsBW02uPy8kBlv/2EPQoBFFwj0GwZ+5xpEwnaw6KhMpKaVwWQrSc2wGBmcuYea/p+QUAhwF0BZANYJVps1UA7nX3WM7kdPkUhUiFsfRMUI9BUJOTA8TEAIDtROaSzkII0VIebUMgohQAAwHsAtB4zjihAAAXLUlEQVSZmctMq04D6OzkPVOIqICICsrLy1t8TIMBSHluojIGYUjnkLwqlsZlIYQneCwgEFEMgA8B/JaZbWbnNs3oo9q/lZmXM3MmM2d27NixRce0jEGojlfGIJyikKwqkcZlIYQneCQgEFEElGBgYOa1psVniKiLaX0XAD964ljWgn0eBFepprOIknQWQoiW8UQvIwLwFoDDzLzQatXHACabnk8G8JG7x7InVSUKx3QWFYhmmUBHCNEynrhDuAnAwwBuI6J9psddABYAGElExwDcbnrtUVJVcoXtBDpXobLuKklnIYRokYBOXWEwAFMerbeZMEanC51up2qcpbOIjwcqKnxeHCGED0jqCpiqSoaugD7sBxAhpMYgOOOsuqyyUu4ShBBNC+iAAAA5tW+h8LbHYDQi5MYgqGmquizUGtuFEC0T2AGBGTh6FOjd298l0YymehYVFcldghDCuYANCAYDkJJsRNhPZ5FiyJMTnUlOjjmVhTppYBZCOBOQAcEyIK04XBmQdvYqOdFZWbxYLemdIhTHaQghXBOQAUEGpDXNPC7BmVAbpyGEcE1ABgQZkNa8nByl15WaUBynIYRoXkAGBBmQ5hr1+RIg8yUIIVQFZEBQzd0TIpPitIT6fAmQ+RKEEKoCMiBYcve0KQHBKAPSmiDzJQghXBWQAQEAciYyCiN7wvi7p2VAWjOkzUUI4YqADQgoL1cuc1NT/V0SzZM2FyGEKwI3IBQWKv+mpPizFAFBrc2FwCgqkvTYQogrAjIgGAxAyt3pyrSZv7pDTmjNsJ4vAWAQjGBTm4KkxxZCmAVcQLCMUq5op4xSLouUE5oLzPMl6PUEtvtvlwZmIQQQgAFBRim7RxqYhRDOBFxAkBOae6SBWQjhTMAFBDmhuUd99DLL6GUhROAFhLvuUiaStyajlF2nPnqZZPSyECKwAoLBAKxaBTBfGXVLBEyeLAPTWkJGLwsh1ARUQFBrUGYG1q/3T3kCmbTFCCHsBVRAkJOY5zhrc2GWwWpChKqACgjSoOw5zlJjAzJYTYhQFVABIS8PiIy0XRYZKQ3KrWE7etmRtCcIEXq8HhCI6E4iOkpEx4noj+7uj7np18J15tHLROrrpSpOiNDi1YBAROEA/gZgNIC+AB4ior6t3V9uLlBfb7usvl6uZN0l7QlCCMD7dwiDARxn5hPMfBnAGgDZrd2ZNCp7h7QnCCEA7weErgB+sHpdbFrWKtKo7B3SniCEADTQqExEU4iogIgKysvLm9w2Lw/QtW20WSajlD1D2hOECDwGg3lKmEGDPLE/bweEEgDdrV53My2zYOblzJzJzJkdO3Zscmc5OcDyyTugRyGIWOZS9gJnd1thYVJtJISWWKYCKPLcPr0dEL4CcC0RpRJRJIAJAD52Z4c5XT5FIaXBWFsvcyl7gbP2hMZGaUsQwlfMV/5hYU46dtTVIXdWvUPmBnd5NSAwcwOA6QA2ATgM4H1mPujWTgsLga5dHQckCI8wtyeEhzuuk7YEIVqv2ZO81XbmK39mU8eOR+thuH0FMHKkMo+8TodTpSp/pG4i1lBH/szMTC4oKGh6o1tuUb6l7dt9U6gQFRamPsaDCDAafV8eIQKZ+SRvfUWv05mqvEdXAYcPK48jR5Dy2iwU1XRy2Ic+7AcU3vAA0KMH0KMHUpY+jaLKGNPaTDAXOGkBbAFm1sxj0KBB7Ex+PrNez0xoZH27cs7Pd7qp8AC9nlkJCY4PvZ7l+xfCxHJuIud/G/pko/rfUtgp2wVRUUxoVN2WyPG4Op15/SBmD5yD27gdUXzANrqGoehiAqZMUdZJG4J35OU5XtGYmccmAPL9i9Bmf+VfVARMeYKBo0eRk/Qf4OBB4MABnDq1Ffbp5gHglLEb8MorQJ8+ykOvR/I1YaoNxfYdPsx/e7m5HmxY9kRU8dTD2R2Cs6tVvd551BbuM1/5NHWnIEQwcuWqnxsbWZ90Wf1vAyeVJzExzEOHsj6m3OW/Idsrf+Wh0zV9Vw6ggD1wDvZ7ELB+OAsIROonJPtbKOEdzr5/QKqORPBRPyEbOT/vJPOKFczTpzPfdBNzTEwT1TtG5sJCZqOxiX06//txKSBZCamAIHcI/tXUXUJzVy5CaIlb9f3mq/527ZSAMH066+MvtOjKvyUn+ZYIqYDQmlso4Tlq378EZhFonF75LzilXPlPncp8ww1NX/UfPcrc2NjMPn1/bgqpgMBsiq66H5VeRnoJBr6Wn+88IEjVndAatatxpzUN5iv/2Fjm4cNZf1VViy58vHnl7ypPBYTAGodw001A27bAp5/6rlDCIiXFeW8GvV7pmSS9joSvGQxKT5tTp5SeOHfdBaxaZdfnP7wOlxojodbTh8AwHj4K9OwJhIU1PWZAo79vItrDzJnu7sfvye1a5ORJZZSe8AtJky38zX6079SpjqN6X1/GDt2lLzW2RTipj6hM1hPQu7eyU9hm/yVCSOVMC5yAUFsLlJWZU/sJP5A02cKXWnvyZ5W7AABo5HCHCxpn2ZLN2X+NRoRUzrTACQjmugq5Q/ArSZMtfEEtn09LTv5qzFf6oXjl76qACAgGA5ByczLC0IiUp++XagkNcDpZUddG9RVCNMH+bmDGDPdO/vYXLOY7gVC98neV5gOC5Urhx2gwwlB0JkrqqjVArT1Bh4u4q24tUpKNzWZ0FKFJLeOn2t1AZaXr+1Q7+f/613In0Cqe6KrkqYdat1MZlKZd9t3tnhz9PetQ7fc+2UIbHH4fT6r02Y+s5/i2Pznt0qzWxdn+9/Xkk/7v9ulvCJVup5KGOXA465aq1yu35yJ0qHXdJGIwq1X7MNS6g9rT6YDJk4H16690MZWuzgpPdTvVfLbT5GT1k4yzOmzhP84alKWhOfjZjwWornbMlKseDJyLjwdiYuTk70uab0NQrat20lVM+JfThubu2rkLFe5zpTtoZaXr/+fx8aT6N754sTQA+5rmA0JODrD8/9VCj0IQWBqINMxZQ3NeuzwYlle7NH2g0DZXu4M6qwJSawBevFi6g2qGJxoiPPVwmsvowAGlBenvf29tm4vwEYe8Lk9+xvlhk1hHl6SxOcCo5gNykgnUlYc0AHsPQqVRGQCwbh1wzz3Al18CQ4f6vmDCLSmJtSg6E+WwXBqbtcPdfEBqpA3Ad0KmURnAlbOGjFIOSKd+dAwGgDQ2a4XaNJCvL2OHgWBKPqBGNHK4wz6IbHsDmquCJAAEFs23IQBQktpFRwOdOvm7JKIVnDU2x13N0q7gB46jglV6BLUwH5AMBAsOmg8IBgOQsuwZhNVUIyWV5KQRgNQamyNQiwtnG2waJ2UEuud5ukeQs3xAr70mPYKCgabbEAIxL7lQZ1NH3Z1RXVmLyovRDttJu4LnqA4Og2NVkDNq1UDyt6dNITEfQm6u462spFgOTDZJxYoIVZccgwEg7QqtpZYjKPe51ieIk2qg0ORWQCCiV4joCBF9Q0T/JKIOVuueJaLjRHSUiO5ozf5l5GvwcjqILf6ibwsSgFypBpoyuRZFLfg7iY+XaiDh/h3CZgD9mDkDwHcAngUAIuoLYAKAdAB3AniNiBy7JjTD6UlD0lYEPNVBbFSDvIongPvuAw4c8E/BNKb1M4RFOZ0hzNngMDn5C7cCAjP/HzM3mF7uBNDN9DwbwBpmrmPmkwCOAxjc0v3n5QG6trb59SVtRXBQnaZwZQRy8voBW7YA/fsDDz4IHDigWh0SClRHBb/u/gxhUhUknPLE6DZTw/S/AUwyPV9qfm56/RaA+528bwqAAgAFycnJNqPv8vOZ42NqGDAyYOT4eBnZGBIqKphzc5ljYjgfE1kXXhMSo5ztRwbHx7duRLB1ini10cYi+MBDI5VdOdFvAXBA5ZFttU0ugH/iSq8llwOC9cM6dUV+vkru9CA9EQgnKipYf1WV6skuPj6wT3SuzBWgXAi5dvJXmycg0L4T0Xo+CwjN7gB4BMCXAHRWy54F8KzV600AftbcvqwDgkyMI5gdT3TOTpbmE6BWrojN5QCYw8Ntr9jVLnbIzZO/5AgKbZoICFAajA8B6Gi3PB3AfgBtAaQCOAEgvLn9WQcEZycCIk9/lULLnF0YOLtrULur9MbJUu2EHx9/pZrH2e9XF1nP8dHVLn8mSRAnXKGVgHAcwA8A9pker1utywXwPYCjAEa7sj+5QxD21K6mnT/Ur7JdvaJ2dnehdvJ3fufS+nI6C3Jy8hfN8VRA0OxIZRmlLMzUZuNSn4Sd4WomTvvpHCMilF43ly9f2cY8ZaN91k/3qZdTRgaL1gqJkcrRUQzlj4cRHy9/HKHKZpRzodJnXq07ZXy861M02k/nWF9vGwwAJQgsX9bg4WDgfIYw6Q4q/E2T6a+v3B1c+aOtqfFjgYSmmE+S1ncN5rEpjhO72151t1QjWjyesknmQWCAY/nl5C/8TZNVRikpyiAce5L4TDTHlYleWhIkwsOBxsbmt1NjPo55H3q9nPiFdwR1lZHkMBKtZV+99NprjiOif/1rlXTcEUBkpO0ynU6547Df1sycAiLcdBMRH688zMd55x0lIDQ0KP9KSgihdZoMCJLDSHiSK0Hi7beBFSvUE7yZtwWunPzVTvgVFcpD8gGJQKXJKiPpYSSEEK4L6iqjnBylu184NQJghIcrryUYCCGE92gyIBgMSkOgMpk3obFReR0qWS6FEMIfNBkQZKY0IYTwPU0GBOllJIQQvqfJgCC9jIQQwvc0GRDy8gBdlO30fzJTmhBCeJcmAwIAREc0QPIYCSGE72gul9GVMQhXho1KHiMhhPA+zd0hSA8jIYTwD80FBOlhJIQQ/qG5gCA9jIQQwj80FxDy8tQnP5EeRkII4V2aCwgAEB0NSA8jIYTwLU31Mqqqss5yqiSblx5GQgjhG5q6QygpkR5GQgjhL5oKCPaTnJtJDyMhhPA+TQUE+ykMzaSHkRBCeJ+mAkLXrtLDSAgh/MUjAYGI/kBETEQJptdEREuI6DgRfUNE17uyn7g4mSlNCCH8xe2AQETdAYwCYF3TPxrAtabHFADLXNlXVZXMlCaEEP7iiTuE/wEwC8rAAbNsAKtZsRNAByLq0tyOpJeREEL4j1sBgYiyAZQw8367VV0B/GD1uti0rEnSy0gIIfyn2YFpRLQFQKLKqlwAz0GpLmo1IpoCpVoJ4eED0djouI30MhJCCO9rNiAw8+1qy4moP4BUAPuJCAC6AdhLRIMBlADobrV5N9Mytf0vB7AcANLSMvnMaSMu1Vy5cZFeRkII4RutrjJi5m+ZuRMzpzBzCpRqoeuZ+TSAjwH8wtTbaCiA88xc5sp+ZaY0IYTwD2/lMloP4C4AxwFcAvCoK28qKgKMRpkpTQgh/IGYufmtfIQok4ECm2V6PVBY6J/yCCFEICCiPcyc6e5+NDVSWY30MBJCCN/QfECQHkZCCOEbmgoIYXalkR5GQgjhO5oKCHo9oA8vBoGh10sPIyGE8CVNBQQwoDoyTQghhNdpagrNoiLACL3l+ZQpynK5SxBCCO/T1B2C0a4HrCS2E0II39FUQFAj3U6FEMI3NB8QpNupEEL4hqYCQhhs64yk26kQQviOpgICkRHmeXYksZ0QQviWpgKCeepMQBLbCSGEr2kqIFiTHkZCCOFbmg0IgPQwEkIIX9J0QJAeRkII4TuaDQjSw0gIIXxLgwFBps4UQgh/0GBAIOlhJIQQfqDBgCA9jIQQwh80GRAA6WEkhBC+ptmAID2MhBDCtzQZEKSHkRBC+J7mAoJMnSmEEP6hqRnTBg0CCgr8XQohhAhNmrtDEEII4R9uBwQi+g0RHSGig0T0stXyZ4noOBEdJaI73D2OEEII73KryoiIbgWQDWAAM9cRUSfT8r4AJgBIB5AEYAsR9WTmRncLLIQQwjvcvUN4EsACZq4DAGb+0bQ8G8AaZq5j5pMAjgMY7OaxhBBCeJG7AaEngCwi2kVE/yWiG0zLuwL4wWq7YtMyB0Q0hYgKiKigvLzczeIIIYRorWarjIhoC4BElVW5pvfHARgK4AYA7xNRWksKwMzLASw3HesCER1tyfv9JAFAhb8L4QIpp2cFQjkDoYyAlNPTenliJ80GBGa+3dk6InoSwFpmZgC7SZkUOQFACYDuVpt2My1rzlFmznRhO78iogIpp+dIOT0nEMoISDk9jYg80mHf3SqjfwG4FQCIqCeASCjR9GMAE4ioLRGlArgWwG43jyWEEMKL3B2YtgLACiI6AOAygMmmu4WDRPQ+gEMAGgBMkx5GQgihbW4FBGa+DGCSk3V5AFqakWi5O+XxISmnZ0k5PScQyghIOT3NI+Uk5YJeCCFEqJPUFUIIIQD4KCAQ0Qoi+tHU1qC2nohoiSnVxTdEdL3VuslEdMz0mOzncuaYyvctEX1BRAOs1hWalu/zVIu/G+UcTkTnTWXZR0QvWq2705RO5DgR/dHP5ZxpVcYDRNRIRHGmdT75PomoOxFtI6JDpvQrM1S28fvv08Vy+v336WI5/f77dLGcWvh9RhHRbiLabyrnn1S2aUtE75m+s11ElGK1rmUphJjZ6w8ANwO4HsABJ+vvArABAEEZ07DLtDwOwAnTv1ebnl/tx3LeaD4+gNHmcppeFwJI0Mj3ORzAOpXl4QC+B5AGpUfYfgB9/VVOu23vAfCpr79PAF0AXG96HgvgO/vvRAu/TxfL6fffp4vl9Pvv05VyauT3SQBiTM8jAOwCMNRum6kAXjc9nwDgPdPzvqbvsC2AVNN3G97U8Xxyh8DM2wFUNbFJNoDVrNgJoAMRdQFwB4DNzFzFzGcBbAZwp7/KycxfmMoBADuhjK/wORe+T2cGAzjOzCdY6RCwBsp37xUtLOdDAP7urbI4w8xlzLzX9PwCgMNwHFXv99+nK+XUwu/Txe/TGZ/9PltRTn/9PpmZq00vI0wP+4bfbACrTM8/ADCCiAitSCGklTYEZ6kuXE6B4QePQ7lqNGMA/0dEe4hoip/KZO1nptvMDUSUblqmye+TiHRQTqQfWi32+fdputUeCOUqzJqmfp9NlNOa33+fzZRTM7/P5r5Pf/8+iSiciPYB+BHKBYjT3yczNwA4DyAerfg+NTVBTqAgJcvr4wCGWS0exswlpGR83UxER0xXyP6wF4CemauJ6C4oAwiv9VNZXHEPgB3MbH034dPvk4hioPzB/5aZf/LWcdzlSjm18Ptsppya+X26+P/u198nK2O4riOiDgD+SUT9mFm1Xc5dWrlDcJbqorUpMLyGiDIAvAkgm5krzcuZucT0748A/gk/Zndl5p/Mt5nMvB5ABBG5k1LE2ybA7nbcl98nEUVAOSkYmHmtyiaa+H26UE5N/D6bK6dWfp+ufJ8mfv19Wh3zHIBtcKyWtHxvRNQGQHsAlWjN9+ntRhGrho8UOG8EvRu2jXa7TcvjAJyE0mB3tel5nB/LmQylHu5Gu+XtAMRaPf8CwJ1+LGcirowxGQzglOm7bQOl4TMVVxrt0v1VTtP69lDaGdr54/s0fS+rASxqYhu//z5dLKfff58ultPvv09XyqmR32dHAB1Mz6MBfAZgjN0202DbqPy+6Xk6bBuVT6CZRmWfVBkR0d+h9CxIIKJiALOhNI6AmV8HsB5KT47jAC4BeNS0roqIXgLwlWlXc9n2ts3X5XwRSt3ca0qbDRpYSXzVGcqtHKD8qN9l5o1+LOf9AJ4kogYANQAmsPILaSCi6QA2QenRsYKZD/qxnAAwDsD/MfNFq7f68vu8CcDDAL411dMCwHNQTq5a+n26Uk4t/D5dKacWfp+ulBPw/++zC4BVRBQOpUbnfWZeR0RzARQw88cA3gLwDhEdhxK8Jpg+Q4tTCMlIZSGEEAC004YghBDCzyQgCCGEACABQQghhIkEBCGEEAAkIAghhDCRgCCEEAKABAQhhBAmEhCEEEIAAP4/mYNGFu+rPm0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "\n",
    "def f(x):\n",
    "    return (100 / x**2) * math.sin(10 / x)\n",
    "\n",
    "fig = plt.figure()\n",
    "\n",
    "xf = np.linspace(1, 3)\n",
    "yf = [f(x) for x in xf]\n",
    "plt.plot(xf, yf, 'r-', label='f=100 / x**2 * sin(10 / x)')\n",
    "plt.legend()\n",
    "plt.xlim(1, 3)\n",
    "\n",
    "text_pt = plt.text(2.05, 30, '', fontsize=12)\n",
    "\n",
    "adaptsim_points, = plt.plot([1], [f(1)], 'bo', label='adaptsim points')\n",
    "\n",
    "def update_adaptsim(i):\n",
    "    eps = 10 ** (-int(i))\n",
    "    \n",
    "    I, m, p = adaptsim(f, 1, 3, eps=eps)\n",
    "    yp = [f(x) for x in p]\n",
    "    \n",
    "    adaptsim_points.set_data(p, yp)\n",
    "    text_pt.set_text(f\"eps={eps}\\npoints={m}\\nI={I: .10f}\")\n",
    "    \n",
    "    return adaptsim_points,\n",
    "    \n",
    "ani = animation.FuncAnimation(fig, update_adaptsim, np.arange(0, 8), interval=200, blit=True)\n",
    "\n",
    "ani.save('adaptsim.gif', writer='imagemagick', fps=2)\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def asr(f, a, b, eps=1e-8):\n",
    "    \"\"\"自适应 Simpson 求积的递归实现\n",
    "    \n",
    "    Args:\n",
    "        f: 要求积的函数\n",
    "        a, b: 求积区间\n",
    "        eps: 目标精度，达到则停止，返回积分值\n",
    "\n",
    "    Returns: \n",
    "        I: 计算结果：积分的近似值\n",
    "    \"\"\"\n",
    "    def simpson(f, a, b):\n",
    "        return (b - a) * (f(a) + 4 * f((a + b) / 2) + f(b)) / 6\n",
    "    \n",
    "    def asrp(f, a, b, eps, sim):\n",
    "        mid = (a + b) / 2\n",
    "        L = simpson(f, a, mid)\n",
    "        R = simpson(f, mid, b)\n",
    "\n",
    "        if (abs(L + R - sim) <= eps):\n",
    "            return L + R + (L + R - sim) / 15\n",
    "        \n",
    "        return asrp(f, a, mid, eps/2, L) + asrp(f, mid, b, eps/2, R)\n",
    "    \n",
    "    return asrp(f, a, b, eps, simpson(f, a, b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7468241328125185"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "asr(lambda x: math.exp(-x ** 2), 0, 1)\n",
    "# actual result (by sympy): sqrt(pi)*erf(1)/2 = 0.746824132812427"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9460830703671792"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "asr(lambda x: math.sin(x) / x, 1e-64, 1)\n",
    "# actual result (by sympy): Si(1) = 0.946083070367183"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.42602475634631"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "asr(lambda x: (100 / x**2) * math.sin(10 / x), 1, 3)"
   ]
  }
 ],
 "metadata": {
  "file_extension": ".py",
  "kernelspec": {
   "display_name": "Python 3.7.6 64-bit",
   "language": "python",
   "name": "python37664bit8b61c5f9b63a430188c2e600c83aebf6"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "mimetype": "text/x-python",
  "name": "python",
  "npconvert_exporter": "python",
  "pygments_lexer": "ipython3",
  "version": 3
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
